Multi-objective optimization of cyclone separators in series based on computational fluid dynamics

Powder Technology 325 (2018) 452–466

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

Multi-objective optimization of cyclone separators in series based on computational fluid dynamics Rafaello D. Luciano, Bárbara L. Silva, Leonardo M. Rosa, Henry F. Meier* Chemical Engineering Department, University of Blumenau, Blumenau 89030-080, Brazil

A R T I C L E

I N F O

Article history: Received 23 July 2017 Received in revised form 30 October 2017 Accepted 13 November 2017 Available online 21 November 2017 Keywords: Cyclones Multi-objective optimization Computational fluid dynamics COBYLA method Gas-solid flow Eulerian-Eulerian modeling

A B S T R A C T The pressure drop and collection efficiency are generally the most relevant parameters when evaluating the performance of gas-solid cyclone separators. In optimization procedures, mathematical cyclone models are used to find the best possible compromise between the pressure drop and the collection efficiency, based on the design preferences. Although several geometry optimization studies on single cyclone configurations have been reported, multiple cyclones in series remain relatively unexplored in this regard. Highly sophisticated models are described in the literature, but researchers mostly use lower fidelity models due to the shorter computational times involved. Also, when more sophisticated models are used, the authors tend to employ surrogate modeling approaches, which reduce the computational times but add another layer of uncertainty. In addition, several authors have optimized either the efficiency or the pressure drop, usually by adding a constraint on the other parameter. In contrast, herein, we describe the multi-objective optimization of three cyclones in series based on high fidelity computational fluid dynamics (CFD) cyclone modeling. A fully automated methodology solving this problem is evaluated, with the COBYLA optimization method and an Eulerian-Eulerian six-phase two-way gas-solid CFD approach simultaneously minimizing the pressure drop and maximizing the efficiency. No surrogate modeling is used here, which means that the results of the CFD simulations are directly used in the optimization procedures. The results obtained indicate that the optimized trios of cyclones outperform the conventional Stairmand (high efficiency) and Lapple (moderate pressure drop) geometries. The minimization of the emission resulted in 30 times less solids being emitted compared with the Stairmand trio. When minimizing the energy consumption, the pressure drop was 50% that of the Lapple geometry. When balancing the pressure drop and emission several results were considerably better than both geometries simultaneously (e.g. 33% less emission than Stairmand combined with 30% less pressure drop than Lapple). This demonstrates that the multi-objective optimization of cyclones in series delivers excellent results and is highly feasible in industrial timescales without compromising the fidelity of the mathematical cyclone model used. © 2017 Elsevier B.V. All rights reserved.

1. Introduction In the class of gas-solid separators, cyclones are the most widely used equipment in industry [1]. These separators mainly rely on centrifugal forces to separate solid particles from gaseous streams. The main advantages of cyclones are that they can be designed for a wide range of operating conditions, and are economical in terms of investment and maintenance costs due to their simplicity of construction. Two conventional geometries are frequently cited when benchmarking new cyclone geometry proposals: the Stairmand [2] and

* Corresponding author. E-mail address: [email protected] (H.F. Meier).

https://doi.org/10.1016/j.powtec.2017.11.043 0032-5910/© 2017 Elsevier B.V. All rights reserved.

Lapple [3] geometries. The former is a high efficiency and high pressure drop design, while the latter sacrifices efficiency for a lower pressure drop. In the literature, cyclone geometrical variables (Fig. 1) are usually shown in the form of dimensionless relations as a function of the barrel diameter (Dc) which is also known simply as the cyclone diameter. The ratios of the Stairmand and Lapple geometries are shown in Table 1. These designs were obtained through physical experimentation and empirical models which predicted the pressure drop and cut-off diameter as a function of the geometrical variables. Cyclones are generally operated in series when the efficiency of a single cyclone is not high enough for the process being studied, with the main drawback being an increased pressure drop and therefore higher energy consumption (in addition to installation costs). The literature on cyclones operating in series is relatively scarce. The main studies focused on this issue are not optimization studies but

R. Luciano et al. / Powder Technology 325 (2018) 452–466

Ds

Lc

Le

Ls

b

Lco

Dc

Dl Fig. 1. Variables of the tangential cyclone separator geometry.

either validation studies [4] or comparison studies [5–8]. The latter are based on physical experimentation tests, and up to four conventional geometries in series, with specific ratios of barrel and cone lengths to body diameter, were evaluated. Some of the conclusions drawn by the authors were that the cyclones in series provided substantially higher collection efficiencies (and shifted the particle size distribution toward smaller sizes) while also increasing the pressure drop, but not linearly (e.g. two cyclones in series did not double the pressure drop). Whitelock and Buser [8] also noted that the use of three or four cyclones instead of two in their process only slightly increased the efficiency, while considerably increasing the static pressure drop. Vegini et al. [4], on the other hand, used numerical CFD simulations with the CYCLO-EE5 code, which is the code used in the present optimization study, and compared the results to experimental data on cyclones in series used in the cement industry. The authors concluded that the code results showed good agreement with experimental data and thus the code has considerable potential for the prediction of the collection efficiency. Among the first studies related to cyclone optimization found in the literature is that carried out by Gerrard and Liddle [9], who optimized the number of cyclones operating in parallel using the conventional Stairmand geometry and empirical modeling. It took a couple of decades until authors started publishing studies on the optimization of the cyclone geometry, with one of the first examples being the work of Ramachandran et al. [10], who aimed to increase the cyclone efficiency while maintaining a restriction on the pressure drop. Later, multi-objective studies aimed at maximizing the efficiency while simultaneously minimizing the pressure drop began to be reported [11]. Currently, in some cyclone optimization studies, the approach was still to select the best out of a few geometries (e.g. Refs. [12,13]). These studies are undoubtedly related to optimization, but in a Table 1 Stairmand and Lapple geometrical ratios. Geometry

Ds/Dc

Dl/Dc

Le/Dc

Ls/Dc

Lc/Dc

Lco/Dc

b/Dc

Stairmand Lapple

0.50 0.50

0.37 0.25

0.50 0.50

0.50 0.62

1.00 1.38

2.50 2.00

0.20 0.25

453

minimalistic sense [14]. In the present study, the optimization is not limited to a simple improvement by choosing the best among three or four cases that were tested. Instead the best possible configuration within existing limitations is sought, usually requiring hundreds of evaluations. The solids collection efficiency and pressure drop are the most commonly evaluated parameters in optimization studies on cyclone performance. The collection efficiency is defined as the proportion of solids that is collected through the underflow in relation to the mass rate of solids that was fed into the cyclone. The pressure drop is the total energy loss of the flow from the inlet to the top outlet. Both of these performance parameters are dependent on the geometrical variables of the cyclone. Ideally, a higher efficiency and lower pressure drop indicate a better geometry. However, in practice, the efficiency and pressure drop parameters are in conflict, and therefore optimization studies on cyclones focus on finding the best possible compromise (or compromises) between these parameters, in accordance to the design priorities. In order to evaluate the numeric values of these parameters in optimization procedures, either accurate mathematical models or a large number of physical experiments are required. In reality, the costs associated with physical experimentation for optimization studies are prohibitive, and therefore authors use sufficiently accurate mathematical models in its place (e.g. Refs. [15–17]). In this context, the empirical and semi-empirical mathematical models are relatively simple and computationally inexpensive. However, they mainly involve the mathematical regression of empirical data, which generally limits their extrapolation capacity. In addition, these models frequently consider only a portion of the variables that constitute the cyclone geometry and operating conditions, and they are limited to the type of geometry for which they were developed, hindering exploration of novel geometries. Nevertheless, due to its ease of use and low computational cost, empirical modeling is widely reported in cyclone optimization literature [10,18,19], and more recently [16,20,21]. There are also phenomenological models which describe the physical phenomena associated with fluid flow, generating a set of equations that can be solved by numerical approaches, such as the computational fluid dynamics (CFD) technique. CFD has the main advantage of accurately predicting important fluid flow characteristics and cyclone performance variables [22], without being restricted by empirical data ranges. This technique also considers all geometrical and operating variables, and can be used to evaluate unique and/or novel geometries (e.g. without a conical region). Its main drawback is the high computational costs compared to empirical modeling, which limits the extent to which it can be applied in optimization studies [14]. Several authors have used CFD techniques in cyclone geometry optimization. One of the earliest examples was described by Hoekstra et al. [23], who applied CFD techniques to obtain data which was used to build a polynomial meta-model. Meta-models, also known as surrogate models, are simplified models of actual models: CFD alone is a model or an abstraction of a phenomenon in the real word, and a meta-model adds another layer of abstraction. Various recent investigations have used meta-modeling approaches in order to optimize the cyclone geometry. These approaches generally fit simpler models, such as polynomials, to the CFD or experimental data through response surface methodology [24–26] or artificial neural networks [21,24,27] techniques. Then, the optimization is performed with the simpler model effectively saving time and computational resources. The main drawback of the meta-modeling approach is that another error source is introduced into the calculations. Safikhani et al. [27], for example, encountered an error of up to 8.8% in their optimized result when comparing the previously regressed CFD meta-model with new CFD

454

R. Luciano et al. / Powder Technology 325 (2018) 452–466

simulations. It should be noted that the authors did not evaluate the uncertainty due to discretization, which is another source of error that could compound that of the meta-model. Regarding the complexity of multiphase CFD modeling, authors generally prefer simpler models that consider only the effect of the gas phase on the solid phase, disregarding the effect of the presence of particles in the gas phase [12,25,27,28]. This modeling approach is known as one-way gas-solid coupling, and is limited to physical situations in which the solids loading is low enough to neglect its effect (generally referred to as “diluted” flow). Sgrott Jr. et al. [29], on the other hand, were the first authors to use a multiphase two-way gas-solid model in cyclone optimization procedures. In this modeling approach, the solid phases affect the gas phase and vice-versa. Thus, it is not limited to the assumption of very low solids loading or very diluted flow to the same extent as in one-way coupling approaches. Developing the discussion of what is considered “diluted” flow for one-way coupling versus that of two-way coupling, Elghobashi [30] generally recommends one-way coupling for volume fractions 100 to 10,000 times lower than that of two-way coupling. The volume fractions suggested for one-way coupling average around 1 × 10 −7 while for two-way coupling they range from 1 × 10 −5 to around 1 × 10 −3 . For the conditions in the present study, the average particle volume fraction reached as high as 7×10 −4 . This result indicates that even two-way coupling approaches the limit of its applicability in this case, even more so considering that the volume fraction in some regions will be higher than the average. Thus, one-way coupling would most likely constitute an oversimplification. The choice of optimization method has reduced importance in studies that either use meta-models of CFD results or empirical modeling, because these models inherently have very low computational costs. This means the optimization step will not take long to complete regardless of the optimization method used. Therefore, these models have been frequently paired with optimization methods that require a large number of function evaluations (which in turn also provide more detailed results), such as the metaheuristic genetic algorithms [15,17,25,27,31,32]. When coupling optimization methods with CFD solvers, considering that each individual simulation can take several hours to complete, the whole procedure can take much longer depending on the efficiency of the optimization method employed. Sgrott Jr. et al. [29] employed a modified version of the NelderMead method [33] with lower and upper bounds for the variables, as proposed by Box [34], obtaining an optimized geometry for maximum efficiency with a restriction on the pressure drop. The procedure took 15 days of computational time to complete, using the dedicated in-house CFD solver CYCLO-EE5 . Amirante and Tamburrano [28] used the commercial software modeFRONTIER to manage the optimization with the commercial solver ANSYS Fluent combined with a genetic algorithm. The procedure took 40 days to complete and they obtained a full Pareto optimal frontier as a result. Mariani et al. [35] used a hybrid optimization algorithm of the add-on Optimate+ from the commercial STAR-CCM+ CD-ADAPCO CFD package. In contrast to the two previously mentioned studies, in this case only the length and the conicity of the vortex finder were optimized, but the heat-exchange was considered in addition to collection efficiency. A geometry with two cyclones in parallel and another downstream was used. The pressure drop was disregarded during the optimization procedure and only verified for the final geometries. The time required for the procedure was not reported. Two main gaps can be noted in the cyclone CFD-based geometry optimization literature: (I) no reports on studies using the configuration of cyclones in series could be found up to 2016, and although a recent study addressed a hybrid parallel-serial configuration [35], the analysis was limited to the vortex-finder and the minimization

of the pressure drop was not among the objectives; and (II) authors mostly use lower fidelity CFD modeling. Only one article describes the use of a multiphase model in which the particles interfere with the gas flow (two-way coupling). The purpose was to obtain a single geometry for maximum efficiency with a restriction on the pressure drop, and the integration of the optimization method with the CFD simulations was manual [29]. In the study reported herein, a six-phase gas-solid two-way CFDbased optimization procedure was applied to the geometries of multiple cyclones associated in series. We used the direct search optimization method proposed by Powell [36], which was found to be more efficient in this case than the approaches based on the Nelder-Mead method described in the literature [20,24,29]. In the optimization procedures, no meta-model nor surrogate model was used. CFD simulations were automatically carried out whenever objective function values were required. The purpose was to evaluate the effectiveness and the methodology of automated cyclone optimization, initially with the multiobjective optimization of one cyclone, followed by the monoobjective optimization of multiple cyclones in series and lastly the multi-objective optimization of multiple cyclones in series. 2. Materials and methods 2.1. Computational fluid dynamics modeling The simulation code CYCLO-EE5 used in this study was developed in-house and employs the finite volume method to solve an EulerianEulerian six-phase model. This code has been the focus of several verification and validation studies (e.g. Refs. [4,22,37-39]), and as a result it has sufficient accuracy for application in optimization studies, as demonstrated by Sgrott Jr. et al. [29]. This code internally uses dimensionless cylindrical coordinates and dimensionless variables with staggered structured grids or meshes in order to calculate the cyclone performance. The staggered meshes are nondimensionalized and generated automatically according to the cyclone geometry, which results in them always having good orthogonal quality. Fig. 2 shows the center of the volumes in a cut of the nondimensionalized pressure mesh for a cyclone minimizing pressure drop. The highlight in the figure also shows the contours of the pressure mesh volumes, in addition to the center of the axial and radial/tangential velocity meshes. Some boundary conditions are introduced with extra volumes outside the domain, and their implementation is dependent on whether the adjacent volumes have faces or the center coincident with the boundary. More details on the numerical grid/mesh can be found in previous publications [22,37]. The modeling procedure in this study was based on previous studies [4,29,38,39], with the partial differential equations solved by the SIMPLEC pressure-velocity coupling method and the k-e -Prandtl turbulence model. The near-wall treatment method is based on the Launder and Spalding [40] wall functions and on the Alexander [41] free vortex function for the tangential velocity. Further details can be found in the previously cited studies. 2.1.1. Multiphase model The following assumptions and hypotheses were considered: a) five solid phases can represent the behavior of the solids carried in the gas phase, with each solid phase having a distinct particle size, density and volume fraction; b) all phases are continuous, interpenetrate one another and have different fluid dynamic properties at the same point in the time-space domain; c) the conservation equations of mass and momentum apply to all five solid phases and to the gas phase, with the flow being incompressible and isothermal; d) turbulence acts only on the gas phase, and can be modeled by a

R. Luciano et al. / Powder Technology 325 (2018) 452–466

455

Fig. 2. Dimensionless grids. The dots (•) represent the volume centers of the pressure grid, the crosses (×) of the axial velocity grid and the empty circles (◦) of the radial/tangential velocities grid.

hybrid anisotropic Reynolds Averaged approach; e) asymmetric disturbance due to the tangential or volute inlet quickly disappears as the flow enters the cyclone, enabling the assumption of axial symmetry with the 3D axisymmetric model; f) drag forces between all phases are responsible for all phase interactions, with the transfer of momentum at the interface between the phases being predicted by a two-way coupling model; g) solid phases are considered hypothetically fluid, due to the physico-chemical interactions between the particles. The strain tensor acting on the solid phases is neglected, with the solid phases obeying the Eulerian momentum equation without pressure and viscous terms. Considering 0 as a generic conservative property of the phase i with f representing its volume fraction and q its density, a general multiphase Eulerian microscopic equation can be written as: ∂ (qf 0) + ∇ • (qf v0) + f ∇ • J 0 = S0 , ∂t

2.1.2. Turbulence modeling As previously mentioned only the gas phase is considered viscous. Therefore all of the equations in this section refer to the gas phase and the subscript g will be omitted. The shear stress model is based on the Boussinesq approximation where, similarly to the general Newtonian fluid model [42], the effective shear stress tensor (T) is proportional to the strain rate tensor (D): T = −2l ef D.

(2)

The effective viscosity is calculated based on the molecular viscosity (l) and the turbulence viscosity (l t ): l ef = l + l t ,

(3)

(1)

where each term is associated with the following physical interpretation:

and the strain rate tensor is calculated as follows, where the superscript T indicates the transpose operation: D=

∂ ∂t (qf 0)

— accumulation of 0 within the control volume; ∇ • (qfv0) — convective transport of 0 through the control volume surfaces; f∇ • J0 — diffusive transport of 0 through the control volume surfaces; S0 — transformation or generation of 0 within the control volume. Several conservative equations can be obtained from Eq. (1). In this study, we considered the conservation of mass, momentum, turbulence kinetic energy and turbulence dissipation rate as detailed in Table 2. Multiphase simulations with one gas phase and five solid phases were performed, totaling six Eulerian phases. In terms of number of equations, a total of 26 conservative partial differential equations are discretized: 6 for gas and solid phases mass, 3 for gas phase momentum, 15 for solid phases momentum, 1 for turbulence kinetic energy and 1 for turbulence dissipation rate.

 1 ∇v + (∇v) T . 2

(4)

As mentioned above and discussed in Meier and Mori [22], a combination of the standard k-e model of Jones and Launder [43] and the Prandtl mixture length model is adopted to calculate the turbulence viscosity, as proposed by Duggins and Frith [44] and Pericleous [45]. According to these authors, the standard k-e model alone is not suitable for cyclones because it is isotropic, which means it adopts the simplification of the viscosity being the same for all directions. The authors therefore proposed that the viscosity in the axial (z) Table 2 Conservative variables of the mathematical model. Variable

0

J0

S0

Mass Gas phase momentum Solid phase i momentum

1 vg vs i

0 T + pI 0

0    qg fg g + ni=1 F drag g,si   qsi fsi g− F drag si,g

Turbulence kinetic energy

k

Turbulence dissipation rate

e

l ef sk l ef sk

∇k

f(G − qe)

∇e

f (C1 G − C2 qe) ke

456

R. Luciano et al. / Powder Technology 325 (2018) 452–466

and radial (r) directions should be calculated through the standard k-e model, while for the tangential or swirl direction (h) the Prandtl model should be used. Considering that vr , vh and vz represent the scalar components of the velocity vector v and l is the mixing length, the turbulence viscosity in the tangential or swirl direction was calculated as: t t t t lh,h = lr,h = lz,h = lin +   ⎫ 12 ⎧    2 z r 2 ⎬ ⎨ ∂v + vrr 2 + ∂v + ∂r ∂z 2 2       . qr l   v ∂v ⎩ 1 r ∂ h 2 + ∂vr + ∂vz 2 + 2 ⎭ h ∂r r ∂z ∂r ∂z 2

where Dh is the hydraulic diameter of the inlet section, and kin is the turbulence kinetic energy at the inlet of the cyclone. The viscosity components associated to the axial and radial directions were calculated according to the standard k-e model: k2 . e

(7)

Representing the turbulence generation and/or dissipation, the turbulence source term G, which appears in Table 2 for k and e, can be calculated in vector notation: G = −T : ∇v,

(8)

or in nondimensional cylindrical notation:      vr 2 ∂vr 2 ∂vz + + Ld 2 ∂r r ∂z   ⎧ ⎫  vh  2 ∂v ∂ + L2d ∂zh 2 + ⎬ 1 ⎨ r ∂r r   + t . ⎭ Rer,r ⎩ Ld ∂vr + ∂vz 2 ∂z ∂r

2 Retr,r



(9)

In order to convert to this form, dimensionless groups are used [22]: Ld =

rref , zref

Retr,r =

(10)

qvref rref l tr,r

,

2.1.3. Multigrid initial condition acceleration Computational fluid dynamics simulations are usually less prone to divergence and can converge faster to a solution depending on the suitability of the initial condition. In this context, a multigrid acceleration technique was employed in order to provide a better initial condition for the fine mesh simulation step. This is similar to strategies reported in the literature by, for instance, Zhang [47] who showed that multigrid methods are efficient for convergence acceleration. In general terms, these methods consist of solving the partial differential equations initially in coarser meshes and then interpolating this solution into finer meshes, considerably reducing the number of internal iterations necessary for convergence. In the case of Zhang [47], the method is inserted in the internal iterative cycles of the code, and not limited to the initial condition. In the present work, a multigrid acceleration technique was employed in order to provide a better initial condition for the simulation with the fine mesh. Fig. 3 shows the effect of improving the initial condition for the Lapple cyclone geometry. It shows the collection efficiency response variable as a function of the number of iterations of the CFD code, for the simulation step in the fine grid. The stability that the better initial condition gives to the fine mesh simulation step is clear: the simulation without the improved initial condition varied from unrealistic and divergence prone −20 % values up to 80%, while with the improved initial condition the simulation was immediately much closer to the actual result and varied by 75.8% to 77.4%. In terms of the number of iterations until convergence, with the improved initial condition the simulation involved 21,357 iterations, while without it 46,926 iterations, that is, more than double the amount, were required. A reduction of 33% in the total simulation time was achieved with the acceleration technique, which involves the gas-only coarse mesh simulation before the solid phases are introduced into the model and the gas-solid coarse and fine mesh iterations. The same comparison was made for the pressure drop and the Stairmand geometry, and the same behavior was observed. Another form of initial condition acceleration used consists of interpolating

(11)

where the reference parameter rref is defined as Dc/2, zref is Lco + Lc + Ls, and vref is the average inlet velocity. The empirical constants used are shown in Table 3, where the value for the mixing length l is based on Dyakowski and Williams [46] and Meier and Mori [22]. The other parameters are at their default values. To summarize, the calculation of the viscosity can be represented by a symmetric tensor analogous to the stress tensor, with the components related to the tangential direction calculated by the Prandtl

80 Collection efficiency (%)

Gadim =

(12)

(5)

(6)

t l tr,z = l tr,r = lz,z = Cl q

⎞ ⎛ t t t t t t ⎞ lr,h l r,z l r,r lke lPr lke ⎜ t t t ⎟ t t t ⎠ ⎝ . ⎝ lh,r lh,h lh,z ⎠ = lPr lPr lPr t t t t t t lke lPr lke l z,r lz,h l z,z ⎛

Further details regarding older versions of the mathematical model can be found in Meier and Mori [22].

t The term lin physically represents an additional turbulence contribution due to inlet effects, which can be calculated according to Pericleous [45]:

t lin = q kin Dh Cl ,

t mixing length theory (lPr , Eq. (5)), and the remaining components t using the standard k-e model (lke , Eq. (7)):

60 40 20 0

With multigrid Without multigrid

− 20 0

Table 3 Constants in the k-e-Prandtl model.

1000

2000

3000

4000

5000

6000

Iteration

Cl

C1

C2

sk

se

l

0.09

1.44

1.92

1.0

1.3

0.034

Fig. 3. Effect of multigrid acceleration on the collection efficiency for the Lapple geometry.

R. Luciano et al. / Powder Technology 325 (2018) 452–466 Table 4 Mesh convergence results. Lapple

Stairmand

DP

g

DP

g

1.4 10.83% 7.70%

2.1 0.79% 0.47%

1.8 5.51% 3.13%

1.9 0.93% 0.51%

the results of a previously simulated cyclone with the most similar geometry to that of the current case, once again in order to provide a better initial condition. This strategy coupled with the multigrid initial condition technique provides even better results in terms of computational costs, contributing to a significant reduction in the time elapsed per simulation to less than 15 min on average for some of the optimization procedures studied. 2.1.4. Grid dependence analysis The grid convergence index (GCI) method, proposed by Roache [48] and further discussed in the literature [49–51], is a method that seeks to provide an estimate of the uncertainty related to the refinement of the numerical mesh in fluid dynamics simulations. It also provides a better understanding of how close the numerical mesh being used is to an ideal numerical mesh with infinite subdivisions or control volumes. In this procedure, the values for the variables of interest are obtained from the simulation with coarse, intermediate and fine meshes. Based on the ratios between the average sizes of the control volumes for each mesh and on the differences between the simulated results, a new value is calculated using the Richardson extrapolation. This value is considered to be an approximation of the ideal value of a numerical mesh with infinite subdivisions or control volumes. Based on this extrapolation, the mesh refinement error is calculated. Finally, with this parameter, the uncertainty estimate is evaluated with the addition of a safety factor. The uncertainty estimated by this method is related to the level of confidence of the results provided by simulations with a given mesh refinement level. Table 4 and Figs. 4 and 5 show the results of the mesh convergence tests performed with the Stairmand and Lapple geometries. Extra mesh refinement points were added in order to confirm that the extrapolation follows the trend of the simulations. Values for the refinement ratios between the coarse and intermediate meshes r32 and intermediate and fine meshes r21 were between 1.3 and 1.4, and are within the method authors’ recommendation of ≥1.3. The result for the apparent or observed convergence order p is on average close to 2, indicating a second-order convergence rate. The

P (Pa)

1280 1240 1200 1160 1120 76.80 (%)

76.45

75.40 0.0000

0.0002 0.0003 1 / Number of volumes

1325 84.8 84.4 84.0 83.6 83.2 82.8 0.0000

0.0001

0.0002 0.0003 0.0004 1 / Number of volumes

0.0005

0.0006

extrapolated values obtained applying the GCI method for the pressure drop seem to overestimate the actual tendency of the redundant points added to the graphs, which indicates that the uncertainty estimated by the GCI method can be considered conservative in this case. Nevertheless, the estimated uncertainty for the pressure drop was higher for the Lapple geometry with approximately 10% uncertainty, while the uncertainty for the collection efficiency was higher for the Stairmand geometry with less than 1% uncertainty. The uncertainty for the intermediate mesh was not much higher than that for the fine mesh, while the computational cost increased disproportionately. The intermediate mesh proportions were therefore used for the optimization procedures, amounting to around 3000 to 5000 three-dimensional axisymmetric control volumes depending on the geometry requested by the optimization procedures. 2.1.5. Convergence criteria A careful evaluation of the convergence in CFD simulations is carried out mainly through assessment of iteration history graphs of residuals and relevant variables. When the user considers that these are sufficiently stable, the simulation is stopped. In practice, this is a subjective criterion and be will dependent on the experience of the individual. In the present study the optimization was fully automated and therefore the first question to be considered was how to instruct the computer when to consider that each simulation is adequately converged. In the area of statistics and probability theories, an appropriate variable to represent what was desired, known as coefficient of variation, is described in Everitt and Skrondal [52]. This coefficient quantitatively measures the variability versus the mean of the results, is limited to positive or absolute scales and noise or oscillations in the results imposes a minimum value that this coefficient will effectively achieve. The limitation to positive or absolute scales is not a drawback here as both the pressure drop and the collection efficiency only have stable physically acceptable values above zero. The coefficient of variation (Cv ) is therefore obtained from the standard deviation (s) and the mean (l) of the results:

Cv =

0.0001

1350

Fig. 5. Mesh convergence study with the Stairmand geometry: • values considered, ◦ redundant points and × extrapolated values.

76.10 75.75

1375

1300

(%)

p GCImedium GCIfine

1400 P (Pa)

Parameter

457

s . l

(13)

0.0004

Fig. 4. Mesh convergence study with the Lapple geometry: • values considered, ◦ redundant points and × extrapolated values.

This coefficient is dimensionless, which represents an advantage (e.g. when compared to the standard deviation) because the same convergence criteria value can be used for both the collection efficiency and the pressure drop, which have different dimensions. The

458

R. Luciano et al. / Powder Technology 325 (2018) 452–466

84.8 Collection efficiency (%)

84.6

Tolerance: 5E-04 Value: 83.74

84.4

2E-04 1E-04 83.68 83.68

84.2 84.0 83.8 83.6 83.4 83.2 0

5000

10000

15000

20000

25000

30000

Iteration Fig. 6. Convergence study in the Stairmand cyclone.

convergence criterion considered in this paper is therefore based on the coefficient of variation. Fig. 6 shows an example of the convergence history of the efficiency for the Stairmand geometry with the coarse grid initial condition (see Section 2.1.3). Variations in the monitored variable are negligible with coefficient values below 2 × 10 −4 , with the same behavior observed for the pressure drop and for the Lapple geometry. Very small coefficient values such as 1 × 10 −3 , are difficult to obtain from the simulations due to the characteristic noise or oscillation associated with the pseudo-transient results. This effectively imposes a minimum value for the coefficient of variation that can be reached without filtering. The simulation is considered stable when this criterion of convergence with the coefficient of variation is satisfied. After this point, the time-averaged values for the relevant variables start being calculated and at intervals of 20 real-time seconds these values are checked (amounting to around 500–1000 iterations). When these time-averages have a relative variation of less than 0.001% between checked points the simulation is stopped. The values for the important variables are then employed in the optimization algorithm. 2.2. Constrained Optimization By Linear Approximation (COBYLA) The COBYLA algorithm was initially proposed by Powell [36] and further detailed in subsequent papers [53,54]. It is a numerical optimization method for constrained problems where the derivative of the objective function is not known or is too computationally expensive to calculate (also known as a derivative-free algorithm). COBYLA incorporates the trust-region concept with consecutive linear interpolations of the objective and constraint functions, while iteratively searching for local minima (or maxima) and reducing the trustregion when no improvements outside of it are found. This method meets the requirements of the current optimization study with the calculation of derivatives being prohibitive and constraints being used. In addition, as will be detailed below, the COBYLA method was found to require considerably less function evaluations for cyclone optimization than the commonly used algorithms based on Nelder and Mead [33]. This is an important point, since CFD simulations with high computational costs are used in this study for the objective function evaluation step. The overall computational time of the optimization procedure will therefore be strongly dependent on the effectiveness of the optimization method used. As previously mentioned, Sgrott Jr. et al. [29] used a variation of the Nelder-Mead method with lower and upper bounds for the variables, in two optimization procedures, one with a 500 mm and the other with a 1000 mm diameter geometry. A hard restriction of 2150 Pa was imposed on the pressure drop and the objective was to maximize the efficiency. Cyclones with a pressure drop higher than

the restriction were discarded. For optimal geometries, in general, the higher the pressure drop, the higher the efficiency will be. In this case, since the pressure drop is disregarded as long as it is below the restriction, it is expected that the optimal result will have a pressure drop very close to that restriction. In Table 5, the optimizations of the previous study are reproduced with the COBYLA method in order to verify the differences in terms of performance and results. In terms of convergence of the optimization method, the advantage of the current study with the COBYLA method is clear. With 30% to 40% less simulations, the result is better than the previous and is considerably closer to the pressure drop restriction with higher efficiencies (i.e. better converged). Due to its excellent performance, the COBYLA optimization method was employed in the present study. 2.3. Desirability function Desirability functions are used when it is desired to operate an optimization study that has multiple important response variables at the same time, with a single objective function. This enables the use of mono-objective optimization methods, which are simpler and more commonly used, in multi-objective studies. To this end, there are unsophisticated methods that are easy to understand and use, and flexible enough to incorporate the preferences of the decision maker. According to a review published by Costa et al. [55], the most popular desirability function method is attributed to Derringer and Suich [56], with extensive practical application. In the area of cyclone optimization, the work of Elsayed and Lacor [25] exemplifies this class of functions in multi-objective optimization. The desirability function method proposed by Derringer and Suich [56] uses predefined parameters that represent how desirable each response variable is, which will be called desirability weights (Wi ). In this way, the values of these weights alter the optimal result that will be obtained by the optimization procedure. According to Fletcher [57], the minimization of a function f is equivalent to the maximization of −f and vice-versa. Using the collection efficiency (g) and pressure drop (DP) as examples, where the objective would be to maximize the g and minimize the DP, the objective function fobj is approximated to a desirability function fdes , with a minimization objective:   fobj = fdes = − (dDP )WDP (dg )Wg 1/(WDP +Wg ) .

(14)

The desirability function therefore represents an estimate of how desirable an evaluation of the results is, by combining the individual desirability (di ) of each response variable and the weights associated. Each di is defined such that 0 ≤ di ≤ 1, with fdes therefore being within the interval [0, 1]. The desirability function also has the property that if any individual desirability is zero (in other

Table 5 Comparison with Sgrott Jr. et al. [29]. Dc = 1000 mm

No. of sim. Ds (m) Dl (m) Le (m) Ls (m) Lc (m) Lco (m) b (m) DP (Pa) g (%)

Dc = 500 mm

Sgrott Jr.

Present

Sgrott Jr.

Present

100 0.506 0.375 0.401 0.711 1.544 2.498 0.207 1941 88.65

63 0.4831 0.390 0.400 0.740 1.660 3.000 0.200 2145 92.37

100 0.246 0.183 0.201 0.372 0.702 1.471 0.103 2075 96.8

72 0.244 0.141 0.200 0.370 0.830 1.500 0.100 2149 97.07

DP > DP max DPmax −DP

DP min ≤ DP ≤ DP max ,

DPmax −DPmin ⎪ ⎪ ⎩1

(15)

DP < DP min

and to the desirability associated with the collection efficiency by ˙ e ): means of the solids mass emission rate (m ˙e>m ˙ e,max m ˙e≤m ˙ e,max . ˙ e,min ≤ m m

˙e ˙ e,max −m m

˙ e,max −m ˙ e,min m ⎪ ⎪ ⎩1

1 µm

0.15 0.10 0.05 0.00 0

DP min = 500 Pa ˙ e,max = rcin m (17)

The variable rcin represents the solids loading ratio at the inlet of the cyclone being optimized, and the emission rates are defined omitting the gas flow rate, since the gas volumetric flow rate (Q) is ˙ = rc × Q ≡ constant throughout the cyclones in series (in reality, m [mass of solids]). The values for the limits were chosen as initial estimates based on previous knowledge and in practice they proved to be surprisingly accurate, since the cyclones simulated remained within these limits. For the optimization of multiple cyclones applying the simultaneous approach, considering Nc as the number of cyclones in series, the following limits were adopted: DP max = (2500 Pa) × Nc DP min = (500 Pa) × Nc   ˙ e,max = (1 − 0.8)(1 − 0.5)(Nc −1) × rcin m

15

20

properties and calculated volume fraction field. The discrete solid phases will change from one cyclone to the next due to the variable collection efficiency depending on the particle size (i.e. larger particles are collected more efficiently than finer ones, skewing the distribution that enters the next cyclone). Several model distribution functions are available and some of them represent the size distributions of several particulate materials very well [1]. Among the most commonly used functions is the log-normal distribution, along the Rosin-Rammler distribution, with applicability to a wide variety of powders. Thus, as commonly reported in the cyclone optimization literature [11,18,29,58], the log-normal distribution model was chosen to represent and constitute the particle size distributions. The PSD for the first cyclone is specified, while for the following cyclones it is calculated based on the results of the previous ones. The solids feed in the first cyclone was considered to have a distribution with a mean diameter of 8 l m and standard deviation of 2.5 l m, resulting in the distribution shown in Fig. 7. The discretized bars and diameters represent the specifications of the phases for the CYCLO-EE5 code. Fig. 8 exemplifies the variation in the regressed particle size distributions from one cyclone to the next in a trio of cyclones in series. Effectively, the average particle size will reduce as the solid stream advances through the cyclones in series. This example was derived from the sequential optimization of multiple cyclones in series with the aim of maximizing the efficiency or minimizing the emission (Section 3.2).

(18)

The global pressure drop for cyclones in series, which will be used in the discussion, is calculated as the sum of the pressure drop (difference between the inlet and overflow pressures) values for all of the cyclones in the series. The global efficiency of cyclones in series is calculated considering the amount of solids emitted from the third cyclone in relation to the amount of solids fed into the first cyclone Once again the values were established based on well-founded assumptions, with the maximum values being analogous to unwanted values, and the minimum values being the utopian or “almost possible” best values. 2.4. Particle size distribution

)

.

10 Particle diameter (µm)

Fig. 7. Particle distribution of the feed.

DP max = 2500 Pa

3

5

(16)

The solids emission rate was used instead of the collection efficiency to standardize and simplify the programming together with the minimization of the pressure drop (the objective of maximizing the efficiency is equivalent to that of minimizing the emission). For the current study, the following limits were adopted for the optimizations of a single cyclone and the sequential optimization of multiple cyclones:

˙ e,min = 200 mg/m m

Log-normal distr. Discrete diameters Discrete phases

0.20

˙e
˙ e,min = (100% − 99.9%) × rcin . m

459

0.25

0.8

1 µm

dg =

⎧ ⎪ ⎪ ⎨0

0.7

Probability density function (

dDP =

⎧ ⎪ ⎪ ⎨0

Probability density function (

words, if the result of one variable is unacceptable), then the global desirability is also zero (the global product is unacceptable) [56]. The definition of the individual desirability functions was also based on Derringer and Suich [56] and adapted to the pressure drop:

)

R. Luciano et al. / Powder Technology 325 (2018) 452–466

1st cyclone 2nd cyclone 3rd cyclone

0.6 0.5 0.4 0.3 0.2 0.1 0.0 0

Another important point to be considered and automated in the current study was related to the particle size distribution (PSD). CYCLO-EE5 works with five discrete solid phases, each with its own

5

10

15

Particle diameter (µm) Fig. 8. Example of the variation of the PSD at the inlet of the cyclones.

20

460

R. Luciano et al. / Powder Technology 325 (2018) 452–466

The log-normal distribution is obtained from the variance (s 2 ) and the mean (l), with the following parameters m and s: ⎛ ⎜ l m = log ⎜ ⎝ 1+  s=

Gas + Solids emitted

⎞

s2 l2

⎟ ⎟, ⎠

(19)

  s2 log 1 + 2 , l

3

Dc = 1 m

(20)

rc = 1 kg/m 3 Q = 6075 m 3/h x50 8 µm

resulting in the probability density function (PDF) as a function of a generic variable x:

1

2

PDF(x) =

log (x)−m 1 −( e 2 ) 2 s . √ e xs 2p

2.5. Operating conditions and variable restrictions For all the CFD-based optimization results presented here, mostly the same operating conditions were used. Exceptions are the number of cyclones in series and which objectives were being evaluated. Fig. 9 shows the general operating conditions that were used. Table 6 shows the restrictions imposed on the seven geometrical variables. The constraints were chosen based on the Lapple and Stairmand geometries with the diameter of the cylindrical section at 1.0 m, as follows:

Dl

Le

Dc = 1 m

(21)

The five discrete phases were regressed in a way that highlights the small particle diameter fractions with more phases (using smaller mass fractions). Smaller particles have lower collection efficiencies and are therefore more relevant for capturing the inefficiency of the cyclones. Thus, representing them with more Eulerian phases is ideal. If the first cyclone was fed with five phases with the same mass fraction we would lose some detail in the distribution of the smaller diameters that are fed into the next cyclone (since the average diameter that represents the smaller phases would be larger due to the larger mass fraction). In this paper, based on previous evaluations a ratio between phase mass fractions of 1.3 was chosen, that is, the solid phase with the smallest mean diameter has a fraction of X%, the next of 1.3 × X% and so on until completing five phases. The emphasis given to the phases with smaller diameters can be observed in Fig. 7. It can also be observed in Fig. 7 that there is an error in the approximation of the continuous log-normal distribution to five discrete phases. In practice, better representations can only be obtained by increasing the number of solid phases in this modeling approach. However, increasing the number of Eulerian phases would bring a series of numerical difficulties related to the very low volume fractions of too many solid phases, along with an increase in the computational cost due to the extra partial differential equations that would have to be solved. Meier et al. [38] demonstrated that a model with three solid phases is more appropriate and accurate than a model with one solid phase for cyclone simulation. Costa et al. [39] further showed that with five discrete solid phases, a good balance is achieved in terms of representation capacity versus computational costs.

Ds

2

Dc = 1 m

— lower limit 20% below those of the Stairmand and Lapple values, and upper limit 20% above; — lower limit 20% below the Lapple value, and upper limit equal to the lower limit of Ds, subtracting 5% of the diameter; — lower limit 20% lower than Stairmand, upper limit equal to the Stairmand value;

Collected solids Fig. 9. Operating conditions.

Ls Lc Lco b

— lower limit equal to the upper limit of Le, upper limit 20% higher than Lapple; — lower limit 20% lower than Stairmand, upper limit 20% higher than Lapple; — lower limit 20% lower than Lapple, upper limit 20% higher than Stairmand; — lower limit 20% lower than Stairmand, upper limit 20% higher than Lapple.

2.6. Optimization process methodology In order for the CFD-based optimization to be fully automated, the CFD simulation procedures also need to be automated and integrated into the optimization method. Fig. 10 shows a simplified flowchart of the main steps necessary in order to automate the CFD-based optimization coupling for singlecyclone optimization. In reality, it was also necessary to deal with diverse situations, e.g., power loss, repeated geometries within the same optimization procedure, failure to write result files, executable exit signal handling, individual monitoring of the simulations and dealing with multiple cyclones in series. The error block indicates the situation when a result was unachievable even after the three different strategy attempts, as shown in the flowchart. To deal with a divergence, the strategy adopted in the algorithm was simply to return a bad objective function value (↑DP and ↓ g) for this geometry to the optimization, in such a way that the optimization method would avoid the region of the parameters close to the geometry that diverged. Initially, the multi-objective optimization of the first cyclone in Fig. 9 was performed, disregarding the existence of the other

Table 6 Restrictions on the geometrical variables. Variable

Lower limit (mm)

Upper limit (mm)

Ds Dl Le Ls Lc Lco b

400 200 400 500 800 1600 160

600 350 500 744 1656 3000 300

R. Luciano et al. / Powder Technology 325 (2018) 452–466

461

Fig. 10. Flowchart representing the simulation-optimization workflow.

The effectiveness of the automated cyclone optimization methodology when approximating the optimal Pareto frontier was first explored for the optimization of a single cyclone (the first one of the series in Fig. 9). In addition, the effects of the geometrical variables on the results and the differences between geometries designed for different purposes were also evaluated. The objective function of this

⎧ Dsinf ⎪ ⎪ ⎪ ⎪ ⎪ Dlinf ⎪ ⎪ ⎪ ⎪ ⎨ Leinf subject to Lsinf restrictions ⎪ ⎪ ⎪ Lcinf ⎪ ⎪ ⎪ ⎪ Lcoinf ⎪ ⎪ ⎩ binf

≤ ≤ ≤ ≤ ≤ ≤ ≤

Ds Dl Le Ls Lc Lco b

≤ ≤ ≤ ≤ ≤ ≤ ≤

(22)

Dssup Dlsup Lesup Lssup , Lcsup Lcosup bsup

(23)

with the functions dDP and dg defined by Eqs. (15) and (16), respectively. The results of the procedures are shown in Fig. 11 and are also compared to those of the Stairmand and Lapple geometries. Ten different weights for the desirability function were used, and they are shown along with the results for the optimizations in Table 7.

100 90 Stairmand

80

Lapple

70 60

Lapple

3.1. Multi-objective optimization of one cyclone

  minimize fobj = − (dDP )WDP (dg )Wg 1/(WDP +Wg ) ,

50 0

500

1000

Stairmand

3. Results and discussion

optimization thus included the minimization of the pressure drop and the emission. The initial geometrical values for all these procedures were the average between the lower and upper bounds. The optimization problem is summarized by:

Collection efficiency (%)

cyclones. The objective function sought to minimize the solids emission and minimize the pressure drop. Two optimization approaches were then evaluated with the mono-objective optimization of multiple cyclones in series, now considering all three cyclones and optimizing the emission but disregarding the pressure drop. One of these was the simultaneous approach, where all cyclones in series are optimized at the same time and considered as a single group of variables. The other is the sequential approach, where each cyclone of the series is optimized alone in the order of the series. These approaches were employed because one of the questions associated with the present multiple cyclones optimization study was whether the optimization procedure would reduce the efficiency in one cyclone to increase the efficiency in the next. This could happen due to the fact that in general, the higher the solids loading is, the higher the efficiency will be, but the constant of this ratio is not known. Therefore, the optimization could predict a reduction in the efficiency, for example, in the first cyclone, in order to increase the solids loading of the second cyclone and gain more efficiency than that which was lost in the first. This would effectively increase the global efficiency of the cyclones in series. To complete the present optimization study, the multi-objective optimization of multiple cyclones in series was investigated, considering all of the information given in Fig. 9. The simultaneous approach was used and the objective function was employed to minimize the global pressure drop and the global solids emission of the three cyclones in series.

Iterations Optimal

1500

2000

2500

Pressure drop (Pa) Fig. 11. Pareto optimal frontier for the first cyclone.

3000

462

R. Luciano et al. / Powder Technology 325 (2018) 452–466

Table 7 Dimensions of the optimized cyclones. Wg

WDP

Ds (m)

Dl (m)

Le (m)

Ls (m)

Lc (m)

Lco (m)

b (m)

g (%)

DP (Pa)

˙ e (g/m3 ) m

0 1 1 2 1 3 2 3 4 1

1 3 2 3 1 2 1 1 1 0

0.600 0.600 0.592 0.574 0.558 0.530 0.518 0.499 0.482 0.400

0.343 0.350 0.342 0.350 0.324 0.349 0.350 0.345 0.349 0.324

0.500 0.400 0.400 0.400 0.400 0.400 0.400 0.400 0.400 0.400

0.744 0.744 0.744 0.715 0.694 0.744 0.698 0.686 0.707 0.664

1.249 1.538 1.351 1.288 1.290 1.310 1.269 1.260 1.295 1.292

2.417 2.778 2.379 2.335 2.335 2.373 2.325 2.330 2.343 2.315

0.300 0.160 0.160 0.160 0.160 0.160 0.160 0.160 0.160 0.160

56.69 88.41 88.04 89.36 90.41 92.26 92.9 93.84 94.78 97.59

591.7 834.3 899.3 992.6 1079 1189 1290 1417 1524 2494

433.1 115.9 119.6 106.4 95.9 77.4 71.0 61.6 52.2 24.1

It can be noted that the procedures resulted in geometries with a better quality than the Stairmand and Lapple geometries in relation to pressure drop and collection efficiency. In fact, four of the ten optimal samples are in the region with a simultaneously higher collection efficiency than Stairmand and lower pressure drop than Lapple. Focusing on the minimization of the emission in relation to the 1 kg/m3 feed, the Stairmand geometry emitted 163 g/m3 while the minimum emission geometry resulted in 24 g/m3 . In terms of collection efficiency and pressure drop, the Stairmand geometry had approximately 83.7% efficiency with a pressure drop of 1330 Pa, while the maximum efficiency obtained was of 97.6% but at a cost of 2494 Pa. Another point worth noting is the dominated result with Wg /WDP = 1/2, while all optimizations ideally should provide non-dominated optimal results. In other words, the result of every optimization should be better in relation to at least one objective when compared to all of the others. However, the resulting pressure drop and collection efficiency are poorer compared to the geometry with Wg /WDP = 1/3. The optimization may have stopped at this non-ideal result due to the limitations of local optimization methods, for instance, these methods are susceptible to stopping in critical points in the objective functions, which include local minima and/or inflection points. Another explanation is the fact that some variables (specifically Ds, Le, Ls and b) are restricted by their bounds, therefore limiting the progress of the optimization. Table 8 shows the correlation coefficients obtained for the variables, considering all iterative and final results for all the optimizations. Further details regarding the Pearson correlation coefficient [59] can be found in Lewicki and Hill [60]. A clear trend can be observed in Table 7 between the optimal diameters of the vortex finder (Ds) and the desirability weights that resulted in these optima. As more emphasis was given to the collection efficiency, the lower the optimal result for Ds became, without exception. In Table 8, this finding is reinforced by the strong correlations between Ds and both the pressure drop and the collection efficiency. The negative sign indicates that an increase in Ds leads to a reduction in the pressure drop and a reduction in the collection efficiency. The absolute value for the coefficient indicates that the relation with the pressure drop shows greater linearity compared to that with the efficiency. These values are in agreement with the behavior observed in Table 7 where initially the decrease in the diameter Ds is gradual as more importance is given to the efficiency when still taking into account the pressure drop. The decrease is accentuated in the last result when the only concern is the efficiency and the pressure drop is disregarded, resulting in a high value for

the pressure drop. This indicates that when the goal is to balance the pressure drop and the efficiency of the cyclone, Ds should not be decreased too much, otherwise the pressure drop will increase disproportionately. The variables of the inlet section, specifically the height (Le) and width (b) of the inlet, also show a more linear correlation with the pressure drop and the efficiency. An increase in either one of those variables tends to result in a decrease in the pressure drop and the efficiency. Both are also more linearly related to the efficiency than to pressure drop. As seen in Table 7, the variables Le and b tended toward their lower bounds in all optimizations, except when the only objective was to minimize the pressure drop. In other words, when the collection efficiency is of some importance to the process, the results indicate that it is generally best to work with the smallest inlet area possible in order to increase the velocity in this section. It is noteworthy however, that erosion effects were not evaluated in this study and these are typically accentuated at high velocities. The values for the geometrical variable Ls were closer to the upper limit when more importance was given to the pressure drop, and they tended to be slightly lower when more weight was given to the collection efficiency. The optimization in which the only objective was to maximize the efficiency, disregarding pressure drop, provided the lowest value (which is still close to the upper bound). This behavior is explained by the correlation coefficients, which indicate a negative linear correlation with the pressure drop, and a weaker negative correlation with the collection efficiency. In other words, increasing the Ls value tends to decrease the pressure drop while also decreasing the collection efficiency, although not as linearly in the latter case. This results in very high values for Ls when balancing the pressure drop and collection efficiency, and lower values when disregarding the pressure drop and maximizing the collection efficiency. The remaining variables did not exhibit such clear trends. Nevertheless, it is also possible to note in Table 7 that all of the optimal values for Dl are close to its upper bound. The correlation coefficients for Dl, shown in Table 8, show a weak negative linear correlation with the pressure drop and almost no correlation with the collection efficiency. This indicates that higher values for Dl tend to reduce the pressure drop while having little effect on the efficiency, which explains its optimal values. In addition, for similar geometrical configurations, Demir et al. [61] observed that, in general, increasing the body and conical heights (roughly equivalent to Lc and Lco, respectively) tends to reduce the pressure drop. This behavior was also observed in the present work through the negative linear correlation values obtained.

Table 8 Correlation coefficients for the relation between the geometrical and the response variables. Dependent variables

Ds

Dl

Le

Ls

Lc

Lco

b

Pressure drop Collection efficiency

−0.9522 −0.5654

−0.2115 0.08518

−0.4404 −0.9175

−0.4798 −0.1888

−0.1416 0.2408

−0.3666 −0.1284

−0.4591 −0.948

R. Luciano et al. / Powder Technology 325 (2018) 452–466

Maximize = 97.59% P = 2494 Pa

Minimize P = 56.69% P = 591.7 Pa

463

In this figure, as expected, the formation of a region of depression in the center and just below the vortex finder is observed. According to Hoffmann and Stein [1], in equipment with rotating flow the greater the distance from the axis of rotation, the higher the pressure will be. This behavior can be observed in all cyclone profiles. In terms of differences between the geometries, the pressure gradients are more pronounced when maximizing the efficiency. With higher pressure gradients, higher velocity values are obtained, which result in stronger centrifugal forces, promoting better separation of the solid phase. When minimizing the pressure drop the pressure gradients are an order of magnitude lower than those when maximizing the efficiency, indicating lower velocities, which will result in lower pressure drop and poorer performance in terms of solids collection.

Balanced = 90.41% P = 1079 Pa

3.2. Mono-objective optimization of multiple cyclones in series The code was then evaluated with the mono-objective optimization of three cyclones in series, according to the operating conditions and constraints shown in Section 2.5. In these optimizations, the sole objective function was to minimize the emission of solids, that is, to maximize the global collection efficiency of the trio. The purpose was also to assess different optimization approaches and the effect of different initial estimates on the optimal results. The optimization problem was then defined as follows: Fig. 12. Comparison of single cyclone optimizations with different objectives.

maximize fobj = gglobal , Fig. 12 shows the geometry results for the optimizations with weights Wg /WDP equal to 1/0, 0/1 and 1/1. In general, a longer geometry tended to result when the objective was to minimize the pressure drop, a behavior similar to that reported in the literature [20,27]. Since CFD simulations were used to model every geometry iteration, it is possible to analyze the results in order to gain a better understanding of the variations in the pressure drop and the collection efficiency resulting from the differences in the geometries. The results for the optimal geometries are therefore compared to those of the Stairmand geometry in Fig. 13, with regard to the pressure field, which is the driving force for the flow inside the cyclone.

2000

1500

500

0

− 500

Stairmand

Minimize

P

Maximize

Fig. 13. Comparison of static pressure fields.

Pressure (Pa)

1000

(24)

with the restrictions of Eq. (23) for each cyclone. The optimization of three cyclones in series was performed applying two different approaches. In the first, named “simultaneous”, the geometrical variables of the series of cyclones are considered to form one group of variables for a single optimization. Thus, in the simultaneous optimization for the three cyclones in series there were 21 variables (7 geometrical variables for each of the three cyclones) and one global objective function. Seven optimizations with different initial estimates were repeated using this strategy. The second approach to optimizing the cyclones in series, named “sequential”, consisted of optimizing the cyclones individually and in order (or sequence). Thus, in the sequential optimization for three cyclones in series, three optimizations were carried out sequentially, each with 7 variables and an objective function. One optimization run applied this strategy. The main drawback of the sequential optimizations approach is that it is unable to identify an optimal configuration that does not consist of the cyclones being progressively as optimal as possible. The possible behavior of reducing the efficiency of one cyclone in order to gain more global efficiency in the trio, mentioned above, was not observed. In the progress of the simultaneous optimization number 5, shown in Fig. 14, from iterations 100 to 150, the efficiencies of the first and second cyclones reach their lowest values, while the efficiency of the third achieves its highest value. However, the efficiency gained in the third cyclone is not sufficient to compensate for the loss in the other cyclones. Therefore, from iterations 150 to around 175, the efficiencies of the first and second cyclones increase while the third loses efficiency. From this point onwards, the optimization tends toward values which are close the others, as seen in Table 9. This was the only optimization in which this behavior was observed. Table 9 also shows that the results of applying the sequential approach were slightly better than those of the simultaneous optimizations with different initial estimates. This could be due to the proximity of the variables to the responses (and the objective function). In other words, the objective functions in the sequential

464

R. Luciano et al. / Powder Technology 325 (2018) 452–466

Stairmand

3.3. Multi-objective optimization of multiple cyclones in series Lapple

The full three cyclones in series multi-objective optimization problem were then evaluated considering all of the operating conditions detailed in Section 2.5. The optimization problem can be summarized as:

Second cyclone

70 Stairmand Lapple

40

  minimize fobj = − (dDP )WDP (dg )Wg 1/(WDP +Wg ) ,

Third cyclone

85 75 65

Stairmand Lapple

55 45 0

50

100 150 200 Optimization iteration

250

300

350

Fig. 14. Progression of a mono-objective optimization of multiple cyclones.

approach were more sensitive to the 7 variables of each cyclone, than the sole objective function was to the 21 variables of the simultaneous approach. This is also related to the way the objective function is defined, which as previously mentioned, differs for the simultaneous and sequential approaches. This slight difference in the results between the sequential and simultaneous approaches is, however, well within the extrapolated uncertainty of up to 1% in the collection efficiency of each individual cyclone. In fact, all of the sequential and simultaneous optimizations with the objective of maximizing the efficiency can be considered sufficiently close and within the uncertainty estimated. A relevant difference between the simultaneous and the sequential approaches was found in the number of simulations required per optimization. The sequential approach required 242 simulations to optimize the trio of cyclones, while the simultaneous approach required an average of 666 simulations. Thus, in the case of the simultaneous approach, it is still debatable whether it is worth the computational time in order to account for possible unusual configurations. Table 9 also shows that the proportional efficiency gain of the optimized geometries versus that of the conventional Stairmand and Lapple geometries increases for each separation stage in the series. For example, on comparing simultaneous optimization 1 with the Stairmand trio, a relative efficiency gain of 17% is observed for the first cyclone, 39% for the second cyclone, and 45% for the third Table 9 Efficiency and emission results for mono-objective optimizations of multiple cyclones.

Lapple Stairmand Sequential Simult. 1 Simult. 2 Simult. 3 Simult. 4 Simult. 5 Simult. 6 Simult. 7

g1 (%)

g2 (%)

g3 (%)

gglobal (%)

Emission (g/m3 )

75.98 83.68 97.73 97.64 97.63 97.66 97.55 97.18 97.79 97.24

54.27 59.37 83.91 82.47 81.98 82.54 80.58 81.98 82.36 81.48

47.7 53.02 79.56 77.01 76.5 76.89 75.09 77.64 76.24 77.92

94.25 96.88 99.93 99.9 99.9 99.91 99.88 99.89 99.91 99.89

57.45 31.15 0.7458 0.9523 1.006 0.9427 1.184 1.136 0.9256 1.127

(25)

with the restrictions of Eq. (23) for each cyclone and the functions dDP and dg defined by Eqs. (15) and (16), respectively. More details about the desirability function can be found in Section 2.3. Fig. 15 summarizes the results for the optimizations performed for three cyclones in series, with the optimal results obtained applying the simultaneous and sequential approaches, and the iterative results for the simultaneous approach. The weights applied to the desirability function used in these optimizations were Wg /WDP = 0/1, 1/4, 1/1, 4/1 and 1/0. The results obtained for the Stairmand and Lapple trios are also reported. Once again, the geometries which simultaneously have a higher collection efficiency than the Stairmand geometry and lower pressure drop than the Lapple geometry were found through the optimization procedure. The sequential approach provided slightly better results for the extremes when the sole objective was either to minimize the pressure drop or minimize the emission. At the intermediate points where the goal was to balance the pressure drop and the emission, it is not possible to directly compare the results for the sequential and simultaneous approaches since the definitions of the desirability and objective functions differ. Fig. 16 was obtained by aggregating the optimal results for the multi-objective optimizations of three cyclones in series with those of the isolated cyclone. Based on this figure, for the operating conditions considered, if the pressure drop in a certain industrial case is limited to below 2000 Pa, a single cyclone with 1 m body diameter is more efficient. On the other hand, if an efficiency ≥96% is required, cyclones in series will be the better option. If efficiencies > 97.5% are required, cyclones in series will be the only option. An explanation is required as to why is it that for pressure drops below 2000 Pa and efficiencies below 96%, a single cyclone separation stage is more efficient and imposes a smaller pressure drop than three cyclone stages. This is possibly due to geometrical constraints and the body diameters being fixed at 1 m, which prohibit

100 98

Stairmand Trio

96 Lapple Trio

94 92 90 88 86 84 1000

2000

Stairmand Trio

55

Lapple Trio

85

Global collection efficiency (%)

Collection efficiency (%)

cyclone. Considering that the progression in separation stages indicates lower solids loading and smaller particle diameters, it can be concluded that there are greater margins for the geometry optimizations with these more severe conditions in cyclones.

First cyclone

100 95 90 85 80 75

3000

Simultaneous Sequential 4000

5000

6000

7000

Global pressure drop (Pa) Fig. 15. Pareto optimal frontier for cyclones in series.

8000

R. Luciano et al. / Powder Technology 325 (2018) 452–466

cyclone alone will be insufficient. In terms of the solids emission for the case studied, to obtain values below 40 g/m3 three cyclones in series are recommended and this is the only choice when the emission needs to be below 24 g/m3 . Three cyclones in series with body diameters of 1 m minimized the emission of solids to 1 g/m3 .

100 Collection efficiency (%)

465

95

90

Acknowledgments

85

The authors would like to acknowledge the support of Petróleo Brasileiro S.A. — Petrobras (cooperation term 0050.0070334.11.9) and CNPq (process number 308714/2016-4).

Three cyclones One cyclone

80 0

1000

2000

3000

4000

5000

6000

7000

Global pressure drop (Pa) Fig. 16. Comparison of the optimizations of single cyclone and three cyclones in series.

the geometries of cyclones in series from achieving higher efficiencies at these flow rates. This behavior could also be explained by the fact that one cyclone is sufficient for these performance conditions and using more cyclone bodies becomes excessive and leads to reduced comparative performance when lower pressure drops are required. 4. Conclusions With the acceleration and interpolation techniques employed, considerably faster optimizations were obtained. The CFD-based optimization code for cyclones in series developed takes 1 to 2 days to optimize a single cyclone, and around 1 week for three cyclones in series. These are industrially viable optimization times. As a basis for comparison, in this study over 10,000 CFD simulations were performed. The results reported herein that the most important geometrical variables in terms of the pressure drop and the collection efficiency are the height and the width of the inlet, and the diameter of the vortex finder. Other variables also have an effect, but are not as noteworthy. Greater optimization improvements are observed when working with more than one cyclone in series, indicating that broader margins are available, under conditions with smaller particle sizes and lower solids loading ratio, when compared to the conventional Stairmand and Lapple geometries. For the case studied, a notable quantity of geometries were simultaneously and considerably better than both the Stairmand and Lapple geometries, highlighting the benefits of performing optimization studies on a case-by-case basis. For the optimization of cyclones in series, if time is an important consideration we would recommend the sequential optimization approach, in which each cyclone is optimized individually and in sequence, instead of working with a single optimization that encompasses multiple cyclones in series together (simultaneous approach). This recommendation is based on the absence in this investigation of any case where the simultaneous approach provided a better result than the sequential approach (even though it is possible that this could be observed for another specific case), while the sequential approach is shown to be nearly three times faster. With the methods described, it was possible to establish ranges for the pressure drop and efficiency where either a single cyclone or cyclones in series are more suitable. For the operating conditions studied, if a pressure drop below 2000 Pa is required, a single cyclone is the better option. Otherwise, if collection efficiencies above 96% are required, three cyclones in series are recommended and this will be the only choice if the efficiency needs to be above 97.5%, as one

References [1] A.C. Hoffmann, L.E. Stein, Gas Cyclones and Swirl Tubes, 2nd ed. ed., Springer-Verlag Berlin Heidelberg. 2010, http://www.springer.com/book/ 9783540746942. https://doi.org/10.1007/978-3-540-74696-6. [2] C.J. Stairmand, The design and performance of cyclone separators, Trans. Inst. Chem. Eng. 29 (1951) 356–383. [3] C.E. Lapple, Processes use many collector types, Chem. Eng. 58 (1951) 144–151. [4] A.A. Vegini, H.F. Meier, J.J. Iess, M. Mori, Computational fluid dynamics (CFD) analysis of cyclone separators connected in series, Ind. Eng. Chem. Res. 47 (1) (2008) 192–200. http://pubs.acs.org/doi/abs/10.1021/ie061501h. https:// doi.org/10.1021/ie061501h. [5] M. Gillum, S. Hughs, B. Armijo, Use of secondary cyclones for reducing gin emissions, Trans. ASAE 25 (1982) 0210–0213. https://doi.org/10.13031/2013. 33505. [6] E. Columbus, Series cyclone arrangements to reduce gin emissions, Trans. ASAE 36 (1993) 545–550. https://doi.org/10.13031/2013.28371. [7] D.P. Whitelock, M.D. Buser, Multiple Series Cyclones for Fine Dust, 2005 ASAE Annual International Meeting, 2005, [8] D.P. Whitelock, M.D. Buser, Multiple series cyclones for high particulate matter concentrations, Appl. Eng. Agric. 23 (2007) 131–136. [9] A.M. Gerrard, C.J. Liddle, The optimal choice of multiple cyclones, Powder Technol. 13 (1976) 251–254. http://www.sciencedirect.com/science/article/ pii/0032591076850103. https://doi.org/10.1016/0032-5910(76)85010-3. [10] G. Ramachandran, D. Leith, J. Dirgo, H. Feldman, Cyclone optimization based on a new empirical model for pressure drop, Aerosol Sci. Technol. 15 (1991) 135–148. http://www.tandfonline.com/doi/abs/10.1080/ 02786829108959520. https://doi.org/10.1080/02786829108959520. [11] G. Ravi, S.K. Gupta, M.B. Ray, Multiobjective optimization of cyclone separators using genetic algorithm, Ind. Eng. Chem. Res. 39 (2000) 4272– 4286. http://pubs.acs.org/doi/abs/10.1021/ie990741c. https://doi.org/10.1021/ ie990741c. [12] M. Wasilewski, L.S. Brar, Optimization of the geometry of cyclone separators used in clinker burning process: a case study, Powder Technol. (2017) http:// linkinghub.elsevier.com/retrieve/pii/S0032591017302310. https://doi.org/10. 1016/j.powtec.2017.03.025. [13] V.V. Biryuk, A.A. Gorshkalev, A.B. Tsapkova, A.A. Shimanov, E.V. Blagin, Optimization of geometric parameters of the cyclone apparatus based on its numerical modeling, J. Phys. Conf. Ser. 803 (2017) 012022. http://stacks.iop.org/1742-6596/803/i=1/a=012022?key=crossref. bc0646810c9f5e6613e1a143cb95124c. https://doi.org/10.1088/1742-6596/ 803/1/012022. [14] D. Thénin, G. Janiga (Eds.), Optimization and Computational Fluid Dynamics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, http://link.springer.com/ 10.1007/978-3-540-72153-6. https://doi.org/10.1007/978-3-540-72153-6. [15] K. Elsayed, Optimization of the cyclone separator geometry for minimum pressure drop using co-kriging, Powder Technol. 269 (2015) 409–424. http:// linkinghub.elsevier.com/retrieve/pii/S0032591014008377. https://doi.org/10. 1016/j.powtec.2014.09.038. [16] A. Alves, J. Paiva, R. Salcedo, Cyclone optimization including particle clustering, Powder Technol. 272 (2015) 14–22. http://linkinghub.elsevier. com/retrieve/pii/S0032591014009140. https://doi.org/10.1016/j.powtec.2014. 11.016. [17] H. Safikhani, Modeling and multi-objective Pareto optimization of new cyclone separators using CFD, ANNs and NSGA II algorithm, Adv. Powder Technol. (2016) http://www.sciencedirect.com/science/article/pii/ S0921883116302278. https://doi.org/10.1016/j.apt.2016.08.017. [18] R.L.R. Salcedo, M.G. Candido, Global optimization of reverse-flow gas cyclones: application to small-scale cyclone design, Sep. Sci. Technol. 36 (2001) 2707– 2731. http://www.tandfonline.com/doi/abs/10.1081/SS-100107221. [19] P.K. Swamee, N. Aggarwal, K. Bhobhiya, Optimum design of cyclone separator, AIChE J. 55 (2009) 2279–2283. http://doi.wiley.com/10.1002/aic.11837. https://doi.org/10.1002/aic.11837. [20] K. Elsayed, C. Lacor, Optimization of the cyclone separator geometry for minimum pressure drop using mathematical models and CFD simulations, Chem. Eng. Sci. 65 (2010) 6048–6058. http://linkinghub.elsevier.com/retrieve/ pii/S0009250910005245. https://doi.org/10.1016/j.ces.2010.08.042. [21] K. Elsayed, C. Lacor, Modeling and Pareto optimization of gas cyclone separator performance using RBF type artificial neural networks and genetic algorithms,

466

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

R. Luciano et al. / Powder Technology 325 (2018) 452–466 Powder Technol. 217 (2012) 84–99. http://linkinghub.elsevier.com/retrieve/ pii/S0032591011005584. https://doi.org/10.1016/j.powtec.2011.10.015. H.F. Meier, M. Mori, Anisotropic behavior of the Reynolds stress in gas and gas-solid flows in cyclones, Powder Technol. 101 (1999) 108–119. http://www. sciencedirect.com/science/article/pii/S0032591098001624. https://doi.org/10. 1016/S0032-5910(98)00162-4. A.J. Hoekstra, J.J. Derksen, H.E.A. Van, Den Akker, Cyclone optimisation based on CFD predictions, Proceedings of the 5th International Conference on Cyclone Technologies, 2000, pp. 53–64. K. Elsayed, C. Lacor, Modeling, analysis and optimization of air cyclones using artificial neural network, response surface methodology and CFD simulation approaches, Powder Technol. 212 (2011) 115–133. http://linkinghub.elsevier. com/retrieve/pii/S0032591011002361. https://doi.org/10.1016/j.powtec.2011. 05.002. K. Elsayed, C. Lacor, CFD modeling and multi-objective optimization of cyclone geometry using desirability function, artificial neural networks and genetic algorithms, Appl. Math. Model. 37 (2013) 5680–5704. http://linkinghub. elsevier.com/retrieve/pii/S0307904X12007263. https://doi.org/10.1016/j.apm. 2012.11.010. X. Sun, S. Kim, S.D. Yang, H.S. Kim, J.Y. Yoon, Multi-objective optimization of a Stairmand cyclone separator using response surface methodology and computational fluid dynamics, Powder Technol. (2017) http://www. sciencedirect.com/science/article/pii/S0032591017305338. https://doi.org/10. 1016/j.powtec.2017.06.065. H. Safikhani, A. Hajiloo, M. Ranjbar, Modeling and multi-objective optimization of cyclone separators using CFD and genetic algorithms, Comput. Chem. Eng. 35 (2011) 1064–1071. http://linkinghub.elsevier.com/retrieve/pii/ S0098135410002656. https://doi.org/10.1016/j.compchemeng.2010.07.017. R. Amirante, P. Tamburrano, Tangential inlet cyclone separators with low solid loading: design by means of 3D fluid dynamic optimization, Eng. Computat. 33 (2016) 2090–2116. http://www.emeraldinsight.com/doi/abs/10.1108/ EC-07-2015-0191. https://doi.org/10.1108/EC-07-2015-0191. O.L. Sgrott, Jr, D. Noriler, V.R. Wiggers, H.F. Meier, Cyclone optimization by COMPLEX method and CFD simulation, Powder Technol. 277 (2015) 11– 21. http://linkinghub.elsevier.com/retrieve/pii/S0032591015001606. https:// doi.org/10.1016/j.powtec.2015.02.039. S. Elghobashi, On predicting particle-laden turbulent flows, Appl. Sci. Res. 52 (1994) 309–329. http://link.springer.com/10.1007/BF00936835. https://doi. org/10.1007/BF00936835. S.I. Pishbin, M. Moghiman, Optimization of cyclone separators using genetic algorithm, Int. Rev. Chem. Eng. (I. RE. CH. E.) 2 (2010) http://profdoc.um.ac.ir/ articles/a/1019360.pdf. N. Fathizadeh, A. Mohebbi, S. Soltaninejad, M. Iranmanesh, Design and simulation of high pressure cyclones for a gas city gate station using semiempirical models, genetic algorithm and computational fluid dynamics, J. Nat. Gas Sci. Eng. 26 (2015) 313–329. http://linkinghub.elsevier.com/retrieve/pii/ S1875510015002759. https://doi.org/10.1016/j.jngse.2015.06.022. J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313. http://comjnl.oxfordjournals.org/content/7/4/308. https:// doi.org/10.1093/comjnl/7.4.308. M.J. Box, A new method of constrained optimization and a comparison with other methods, Comput. J. 8 (1965) 42–52. http://comjnl.oxfordjournals.org/ content/8/1/42.short. F. Mariani, F. Risi, C.N. Grimaldi, Separation efficiency and heat exchange optimization in a cyclone, Sep. Purif. Technol. (2017) http://linkinghub.elsevier. com/retrieve/pii/S1383586616323656. https://doi.org/10.1016/j.seppur.2017. 02.024. M.J. Powell, A direct search optimization method that models the objective and constraint functions by linear interpolation, Advances in Optimization and Numerical Analysis, Springer. 1994, pp. 51–67. http://link.springer.com/ chapter/10.1007/978-94-015-8330-5_4. H.F. Meier, J.J.N. Alves, M. Mori, Comparison between staggered and collocated grids in the finite-volume method performance for single and multi-phase flows, Comput. Chem. Eng. 23 (1999) 247–262. http://www. sciencedirect.com/science/article/pii/S0098135498002701. https://doi.org/10. 1016/S0098-1354(98)00270-1. H.F. Meier, A.A. Vegini, M. Mori, Four-phase Eulerian-Eulerian model for prediction of multiphase flow in cyclones, J. Comput. Multiphase Flows 3 (2011) 93–105. http://cmf.sagepub.com/content/3/2/93. https://doi.org/10. 1260/1757-482X.3.2.93. K.K. Costa, O.L. Sgrott, Jr, R.K. Decker, E.L. Reinehr, W.P. Martignoni, H.F. Meier, Effects of phases numbers and solid-solid interactions on the numerical simulations of cyclones, Chem. Eng. Trans. 32 (2013) 1567–1572. http://www. aidic.it/cet/13/32/262.pdf. https://doi.org/10.3303/CET1332262.

[40] B.E. Launder, D.B. Spalding, The numerical computation of turbulent flows, Comput. Methods Appl. Mech. Eng. 3 (1974) 269–289. http://www. sciencedirect.com/science/article/pii/0045782574900292. https://doi.org/10. 1016/0045-7825(74)90029-2. [41] R. Alexander, Fundamentals of cyclone design and operation, Proc. Aus. Inst. Min. Met. 152 (1949) 152–153. [42] D.C. Wilcox, Turbulence Modeling for CFD, 3rd ed. ed., D C W Industries, La Cañada, Calif, 2006. [43] W.P. Jones, B.E. Launder, The prediction of laminarization with a two-equation model of turbulence, Int. J. Heat Mass Transf. 15 (1972) 301–314. http://www. sciencedirect.com/science/article/pii/0017931072900762. https://doi.org/10. 1016/0017-9310(72)90076-2. [44] R. Duggins, P. Frith, Turbulence anisotropy in cyclones, Filtr. Sep. Nov/Dec (1987) 394–397. [45] K.A. Pericleous, Mathematical simulation of hydrocyclones, Appl. Math. Model. 11 (1987) 242–255. http://www.sciencedirect.com/science/article/pii/ 0307904X87901399. https://doi.org/10.1016/0307-904X(87)90139-9. [46] T. Dyakowski, R.A. Williams, Modelling turbulent flow within a smalldiameter hydrocyclone, Chem. Eng. Sci. 48 (1993) 1143–1152. http://www. sciencedirect.com/science/article/pii/000925099381042T. https://doi.org/10. 1016/0009-2509(93)81042-T. [47] J. Zhang, Multigrid Acceleration Techniques and Applications to the Numerical Solution of Partial Differential Equations, Doctor of Philosophy, George Washington University. 1997, http://citeseerx.ist.psu.edu/viewdoc/download? doi=10.1.1.110.1949&rep=rep1&type=pdf. [48] P.J. Roache, Perspective: a method for uniform reporting of grid refinement studies, J. Fluids Eng. 116 (1994) 405–413. https://doi.org/10.1115/1.2910291. [49] P.J. Roache, Quantification of uncertainty in computational fluid dynamics, Annu. Rev. Fluid Mech. 29 (1997) 123–160. http://www.annualreviews.org/ doi/abs/10.1146/annurev.fluid.29.1.123. [50] P.J. Roache, Verification of codes and calculations, AIAA J. 36 (5) (1998) 696– 702. https://doi.org/10.2514/2.457. https://doi.org/10.2514/2.457. [51] I.B. Celik, U. Ghia, P.J. Roache, C.J. Freitas, H. Coleman, P.E. Raad, Procedure for estimation and reporting of uncertainty due to discretization in CFD applications, J. Fluids Eng. 130 (2008) 078001. http://FluidsEngineering. asmedigitalcollection.asme.org/article.aspx?articleid=1434171. https://doi. org/10.1115/1.2960953. [52] B. Everitt, A. Skrondal, The Cambridge Dictionary of Statistics, Cambridge University Press, Cambridge; New York, 2010, http://www.books24x7.com/ marc.asp?bookid=36106. oCLC: 668966755. [53] M.J.D. Powell, Direct search algorithms for optimization calculations, Acta Numer 7 (1998) 287. http://www.journals.cambridge.org/abstract_ S0962492900002841. https://doi.org/10.1017/S0962492900002841. [54] M.J. Powell, A view of algorithms for optimization without derivatives, Mathematics Today-Bulletin of the Institute of Mathematics and its Applications 43 (2007) 170–174. http://citeseerx.ist.psu.edu/viewdoc/ download?doi=10.1.1.591.6481&rep=rep1&type=pdf. [55] N.R. Costa, J. Lourenco, Z.L. Pereira, Desirability function approach: a review and performance evaluation in adverse conditions, Chemom. Intell. Lab. Syst. 107 (2011) 234–244. http://www.sciencedirect.com/science/article/pii/ S0169743911000797. https://doi.org/10.1016/j.chemolab.2011.04.004. [56] G. Derringer, R. Suich, Simultaneous optimization of several response variables, J. Qual. Technol. 12 (1980) 214–219. [57] R. Fletcher, Practical Methods of Optimization, 2nd ed. ed., Wiley, Chichester; New York, 2000. [58] R.I. Kermode, A geometric programming solution to the optimal recovery of solids from gas cyclones, Eng. Costs Prod. Econ. 7 (1982) 63–68. http://www. sciencedirect.com/science/article/pii/0167188X82900106. https://doi.org/10. 1016/0167-188X(82)90010-6. [59] K. Pearson, Note on regression and inheritance in the case of two parents, Proc. R. Soc. Lond. 58 (1895) 240–242. http://www.jstor.org/stable/115794. [60] P. Lewicki, T. Hill, Statistics: Methods and Applications, Statsoft, Inc., Tulsa, OK, 2006, http://www.statsoft.com/Textbook. [61] S. Demir, A. Karadeniz, M. Aksel, Effects of cylindrical and conical heights on pressure and velocity fields in cyclones, Powder Technol. 295 (2016) 209– 217. http://www.sciencedirect.com/science/article/pii/S0032591016301437. https://doi.org/10.1016/j.powtec.2016.03.049.