Coupled Computational Fluid DynamicsâDiscrete Element Method

Article pubs.acs.org/crystal

Coupled Computational Fluid Dynamics−Discrete Element Method Simulations of a Pilot-Scale Batch Crystallizer B. Ashraf Ali,†,‡,∥ M. Börner,§,⊥ M. Peglow,§ G. Janiga,† A. Seidel-Morgenstern,‡,§ and D. Thévenin*,† †

University of Magdeburg “Otto von Guericke”, Institute of Fluid Dynamics and Thermodynamics, 39106 Magdeburg, Germany Max Planck Institute for Dynamics of Complex Technical Systems, 39106 Magdeburg, Germany § University of Magdeburg “Otto von Guericke”, Institute of Process Engineering, 39106 Magdeburg, Germany ‡

S Supporting Information *

ABSTRACT: Computational fluid dynamics (CFD) coupled with the discrete element method (DEM) has been used to investigate numerically crystal dynamics in an existing pilot-scale batch crystallizer. The CFD-DEM combination provides a detailed description of crystal dynamics considering a four-way coupling. In a previous analysis,1 CFD had been coupled with a discrete phase model (DPM) using a simple one-way coupling. The corresponding predictions are then compared with those obtained through four-way coupling considering KH2PO4 crystals in water. From the CFD-DEM simulation, it is possible to investigate quantitatively the driving force controlling crystal growth and the interaction of crystals with reactor walls, baffles, and impellers. This delivers essential data for process improvement. Different seeding procedures were also compared. The seed crystals have been injected either within the complete liquid volume or, as in the experiments, through a funnel. By varying the most important crystallization process parameters, we found optimal conditions for a liquid phase volume in the crystallizer of 24 L, for injection through a funnel above the baffle, and for an initial seed crystal size of 0.5 mm.

1. INTRODUCTION Crystallization is an essential separation and purification process in the process industry.2 Typically, the primary goal is to generate pure particles with a set size and shape. Particle characteristics and process parameters affect the properties as well as the course of subsequent downstream processing steps. In order to obtain the desired product quality, attention should be paid to several important operating conditions, such as temperature, supersaturation, mixing quality, mode of operation, geometry of the crystallizer, or seeding strategies.3 The desired crystal size distribution (CSD) can differ depending on the specific requirements of the producer. In bulk chemical processes, large particles are produced to ensure fast separation times, maximizing the economic efficiency of the process. Conversely, in the pharmaceutical industry, smaller particles are usually required to improve the bioavailability of the active ingredient. In either case, there is a need to target a predefined particle size that, if not met, will result in additional processing costs or even in product losses. The fundamental understanding of crystallization has improved significantly over the past 30 years. Despite these improvements, this single unit operation still poses immense challenges when designing and optimizing new chemical processes. In particular, the potential of numerical simulations for design and optimization of crystallizers still appears to be large.4 This is the central issue considered in this article. Crystal growth rates are crucial for the process outcome. The rates of crystal growth depend on both transport rates and rates of surface integration steps. Regarding the focus of this study, it © 2014 American Chemical Society

is imperative to note that the former rates depend on the specific local velocity differences between the fluid and solid phases present in the crystallizer. These velocity differences have a strong influence on the thickness of the boundary layers around the crystals, which need to be passed by the solute molecules in order to be available for surface integration (see section 4). Experimentally, crystal growth rates can be measured directly by observing the growth of a single seed crystal or through quantifying the growth of a suspension of many crystals.5 Crystal growth kinetics is frequently determined by combining experimental observations with the application of population balances.6 Here, the parameters are estimated indirectly based on evaluating measured desupersaturation curves7 or their initial derivatives.8 Seeding is an effective technique to depress spontaneous nucleation and make crystal growth the dominant process. This process allows for a better target of CSD and yield.9 In this sense, seeding is one of the most critical steps in optimizing crystallization behavior. Experiments show that two different seeding strategies may often lead to entirely different crystal populations.3,10 Seed crystals with the correct initial particle size, mass, and crystal form must be introduced at the right time and at the right position. Seeding can also play an important role in crystal morphology.9,11 Large seeds generate more secondary nuclei in agitated systems than do small seeds, due to Received: July 21, 2014 Revised: November 20, 2014 Published: November 21, 2014 145

dx.doi.org/10.1021/cg501092k | Cryst. Growth Des. 2015, 15, 145−155

Crystal Growth & Design

Article

greater contact probabilities and collision energies.10 Increasing the seed mass and decreasing the seed mean size increases the seed surface area for growth, therefore promoting growth over nucleation.10 Also, seeding should take place at the moment when the solubility concentration is achieved. The excessive buildup of supersaturation prior to seeding should be avoided because it may lead to either spontaneous or seed-induced nucleation. Optimal seeding can change a poorly behaved, inconsistent crystallization process, to one that is repeatable and produces particles with the required particle size.12 For the results discussed in this article, the impact of seeding will be evaluated by determining the effect of seed size and injection position on crystal growth. Since the balance between hydrodynamics and kinetics will control the growth process, an accurate numerical modeling of crystallization requires suitable models for both fluid dynamics in the crystallizer and representing the key phenomena, such as nucleation, growth, and breakage.13 Those models then have to be coupled. Two main approaches have been documented in the scientific literature for that purpose, either (1) considering the particle trajectories in a Lagrangian frame, as done in the present work, or (2) through a statistical description of the crystal population in a Eulerian framework.14 The latter approach has been used in most recent research works concerning crystallization, and is more efficient from the computational point of view. However, collision effects must be suitably modeled through corresponding integral terms,15 and accurate slip velocities can only be obtained with more complex, multifluid formulations. The former Lagrangian approach becomes unacceptably expensive in terms of computing time when considering too many particles. This is why merely 10000 crystals are considered in this article, although the real crystallizer contains many more particles. One major advantage of this approach is that it can be used to directly investigate slip velocities and collision probabilities. The major drawback is its high computational cost, usually limiting the application of CFD-DEM simulations to relatively small-scale installations with not too many particles. Considering a few selected numerical studies of crystallizers, Wei et al.16 analyzed a semibatch precipitator using CFD to predict transient behavior of CSD and supersaturation distribution. In their study, a population balance equation (PBE) was considered and integrated into a CFD solver to describe the generation and transport of the crystal phase. Then, the impact of the feed locations and impeller speed on mixing and CSD in precipitation reactors were investigated. It was found that feeding into the region where mixing is poor may lead to a larger average crystal size with a more widely spread distribution. Liiri et al.17 investigated local hydrodynamics and KDP crystal growth in a 100 dm3 suspension crystallizer equipped with two Rushton turbine impellers. The authors used a combination of CFD with a multiblock model involving 66 compartments and 26 crystal size classes to simulate crystal growth. Compared with experimental measurements, the observed trends were reasonably good, even if differences were observed in some parts of the reactor, perhaps due to the missing model for breakage and nucleation. Rane et al.4,18 recently applied a PBE solver coupled with a CFD code to obtain the CSD for a variety of crystallizers. Eleven industrial crystallizers were considered,4 comparing the width of the obtained CSD and highlighting advantages and drawbacks of the different systems. Using the same approach,

the impact of impeller design and rotation speed was then considered more specifically, in a first effort toward developing corresponding scaling laws.18 To the author’s best knowledge, there are currently no reports in the literature relying on a coupling between computational fluid dynamics (CFD) and discrete element method (DEM) in order to describe crystal dynamics and indirectly crystal growth in a batch crystallizer. With the CFDDEM approach, a four-way coupling is achieved for the twophase flow: hydrodynamics controls to a large extent crystal motion, flow modifications induced by the particles are taken into account, and collisions (crystal−wall and crystal−crystal) are also accounted for. Such a coupled approach may be important since particle dynamics depend directly on mixing and hydrodynamic conditions.16,19,20 At the same time, for high crystal number density (dense conditions), the flow will become influenced by the particles. Crystal collisions, possibly leading to attrition, are also particularly essential for practical considerations. In our earlier work considering the same crystallizer,1 hydrodynamics was investigated in a similar manner using CFD, but crystal dynamics was described by a simple, one-way coupling. However, collisions cannot be accurately modeled with the aforementioned one-way coupling, although they play an essential role for the process. In order to overcome this limitation, the impact of the most important process parameters is analyzed with a coupled CFD-DEM approach in this study. The article is organized as follows. First, the geometry of the considered batch crystallizer is introduced in section 2. In section 3, the applied four-way coupling model is depicted. The connection between hydrodynamics, slip velocity, and crystal growth is discussed in section 4, before presenting and discussing the simulation results concerning volumes, seed sizes, and injection positions in section 5.

2. GEOMETRY OF THE CRYSTALLIZER AND SYSTEM Figure 1 shows the experimental setup used at the Max Planck Institute for Dynamics of Complex Technical Systems in

Figure 1. Geometry of the investigated draft tube crystallizer.

Magdeburg. The total reactor volume is 35 L. Equipped with an impeller, the considered crystallizer is centrally mounted on a shaft. The inner diameter of the crystallizer is 310 mm, while the diameter of the impeller is 180 mm. A stationary baffle (draft tube) begins 90 mm above the bottom of the reactor and extends vertically over 190 mm. The investigations in this article consider mainly a rotational impeller speed of 300 rpm (standard value employed in the experiments). In all 146



Article

need to be fulfilled. Here, ε is the fluid volume fraction, ρ is the fluid density, u is the fluid velocity, p is the pressure shared by all phases, g is the gravitational acceleration, and S is the momentum sink due to interaction with the particles. The second term on the right-hand side of eq 7 corresponds to the momentum flux due to shear stresses and turbulence, the third term represents body forces, and the fourth term is the interaction force between the two phases. In this work, the turbulence of the fluid flow is modeled by the k-epsilon approach. More details concerning CFD can be found in our previous article.1 3.3. CFD-DEM Coupling. A Eulerian−Lagrangian coupling method is used to couple the particle dynamics (DEM) with the fluid flow (CFD). In the coupling between the discrete and the continuous phase, the mutual interactions are accounted for. The interaction force between the two phases is modeled in the continuous phase through a momentum sink term, originating from the drag force that arises because of the slip velocity between the phases. Momentum sink S is calculated by

simulations, KH2PO4 (KDP) crystals are injected in water. The exact reactor geometry has been implemented in the CFD. More details can be found in our previous article.1

3. CFD-DEM COUPLED SIMULATIONS 3.1. DEM Equations and Models. The discrete element method (DEM) is a numerical simulation technique to describe the motion and the interactions of discrete particles in a particular domain, which were first described in the scientific literature at the end of the 1970s. The motion of each particle is calculated applying Newton’s second law of motion. Both translational and rotational motions of each individual particle are determined. In the absence of any continuous flow, they read: mi

Ii

dvi = dt

dωi = dt

∑ (Fijn + Fijt + mig )

(1)

∑ (R i × Fijt − μr R i|Fijn|ωi)

(2)

S=

where indices i, j describe the interacting particles, vi and ωi are the translational and angular velocity of the particles, Ii is their moment of inertia, and Ri is the particle radius. The symbols Fnij and Ftij denote the normal and tangential contact force, with μr being the coefficient of rolling friction. Ii(dωi/dt) is the torque due to tangential stress, and ∑(Ri × Ftij − μrRi|Fnij|ωi) is the rolling friction torque arising from elastic hysteresis. The particle−particle and particle−wall interactions are described by the standard Hertz−Mindlin contact model.21 The contact forces of all interacting partners are a function of the normal contact force, the radii of curvature, and the modulus of elasticity, 4 Fn = E* R *δn 3 3

1 1 1 = + R* Ri Rj

FD = 0.5C DρA|Δv|Δv

(9)

where A is the projected area of the particle and Δv is the relative velocity between the fluid and the particle. The drag coefficient, CD, depends on the local Reynolds number (Re) of the particle and is computed according to eqs 10 and 11: Re =

(3)

αρL|Δv| η

(10)

where η is the fluid viscosity, L is the diameter of the particle’s bounding sphere, and α is the free volume of the CFD mesh cell. The correlation for CD then reads: ⎧ 24 Re ≤ 0.5 ⎪ ⎪ Re CD = ⎨ 0.697 )/Re 0.5 < Re ≤ 1000 ⎪ 24(1.0 + 0.15Re ⎪ ⎩ 0.44 Re > 1000

(4)

(11)

To calculate the solid volume fraction in a CFD mesh cell, a point search algorithm is used where regular sample points are taken within the bounding box of a particle and retained if the point lies within the bounding surface of the particle.22 Each sample point is then checked to determine which CFD mesh cell belongs to the corresponding particle. The particle mass within a particle mesh cell is finally calculated by

(5)

where Ei and Ej are the Young’s moduli of particles i and j, respectively, pi and pj are the corresponding Poisson ratios, and Ri and Rj are the radii of particles i and j. 3.2. CFD Equations and Models. The fluid flow in the system is described by a continuum approach. Hence, the discretized, volume-averaged Navier−Stokes equation is solved. For each finite-volume element of the mesh-based solution, the mass conservation

∂(ερ) + ∇(ερu) = 0 ∂t

(8)

Here, V is the volume of a CFD mesh cell, and ∑FD is the summation of the drag forces exerted on the fluid in the mesh cell. The free stream drag model adopted reads

where Fn is the normal force, E* is the equivalent Young modulus, and δ is the overlap for particles in contact. To predict the equivalent Young modulus, the following equations are used: 1 − pj 2 1 − pi 2 1 = + E* Ei Ej

∑ FD V

m p (t ) =

ρp ncVp

∑ particles

N

(12)

Here, nc is the number of sample points contained within the mesh cell of particle number p, N is the total number of sample points of the particle, ρp is the particle density, and Vp is the volume of the particle. As often done for particulate flows, the lift force acting on the particles is neglected in the present simulation. The possible influence of the lift force is the subject of our current work.

(6)

and the conservation of momentum ∂(ερu) + ∇(ερu ⊗ u) = −∇p + με∇2 u + ερg − S ∂t (7) 147



Article

The CFD is first performed as a single-phase, transient calculation that is iterated until convergence for a specific time step. A drag force is then calculated by taking into account the current position of the DEM particles using the fluid velocity in the mesh cell within which each particle is located. While freezing the flow solution, the DEM software then takes sole control of the simulation and performs iterations on its side. In most cases, the DEM time step is noticeably smaller than the employed CFD time step.23 After the DEM finishes its own iteration process, control is passed back to the CFD. A momentum sink term S is then added to each of the CFD mesh cells containing particles in order to take into account transfer from the fluid to the particles.24 This double iteration process in time is schematically depicted in Figure 2. More details concerning the coupling between CFD and DEM can be found in particular in ref 25.

smaller than CFD time steps, as needed to correctly capture the contact behavior of particles.24 A typical CFD simulation requires roughly 4 days of computing time on a standard desktop PC (16 GB memory, 64 bit, 3.30 GHz) to obtain periodic hydrodynamic conditions without any particle (meaning for the present conditions 5 s of flow time, see ref 1 for more details). The DEM is coupled with the CFD after the fluid flow reaches a periodic steady state in the crystallizer. After that, a fully coupled CFD-DEM simulation is started. The duration of this coupled simulation depends strongly on the number of considered particles (particle loading). For instance, a four-way coupled simulation involving 10000 crystals takes roughly 25 days of computing time on the same PC.

4. CONNECTION BETWEEN HYDRODYNAMICS AND CRYSTAL GROWTH Crystal growth requires both transport and surface integration. First, the solute molecules need to be transported from the bulk of the solution through a boundary layer to the crystal surface by diffusion before they can be attached to the surface. These processes are driven by specific local concentration differences. A more detailed analysis requires knowledge about the specific orders of the aforementioned processes with respect to the driving concentration gradients.2,5,26−29 For the purpose of this study, it is instructive to exploit the simplifying assumption that both processes are linear and that the rates of mass transfer and surface integration are equal. Then, a total (overall) growth rate constant, ktot can be derived:5 k tot =

1 1 kdif

+

1 ksur

with kdif =

D δ

(13)

In the eq 13, kdif is the transport rate constant, ksur is the surface integration rate constant, D is the component diffusivity in the fluid phase, and δ is the boundary layer thickness. This thickness is strongly influenced by the hydrodynamics in the two-phase system, in particular by the local slip velocities.27 Higher slip velocities decrease the layer thickness δ. According to eq 13, this increases kdif and ktot and, thus, favors crystal growth. This trend is generally valid and holds also if more realistic nonlinear (positive) orders are needed to describe the rates of the underlying subprocesses.

Figure 2. Coupled calculation methodology adopted to investigate crystal dynamics.

Crystals do not change size in the present simulations. Consequently, it is not necessary to account for mass exchange between the phases. In the real crystallization process, crystal growth occurs on a time scale of hours, much longer than any CFD-DEM simulation. Therefore, the driving force for crystal growth is determined from the simulations, but there is no direct simulation of the crystal growth itself. Due to isothermal conditions in the crystallizer, energy exchange is neglected in the coupling of both phases. The simulations are performed by using the commercial software packages Ansys-Fluent 12.1 (for CFD) and EDEM 2.3 (for DEM). For both software packages an appropriate coupling module exists. In CFD, a first-order upwind discretization is applied, together with the SIMPLE algorithm for pressure−velocity coupling. Tests at the beginning of this study have shown that a second-order discretization does not lead to significant changes in the solution. A sliding mesh approach is employed to compute the transient behavior of the particle−liquid flow, as described in our previous publication.1 The (fixed) time step used in the CFD computation is 1 ms, while the time step in DEM is 16.1 μs. This is equivalent to 39% of the Rayleigh time, leading to 62 DEM subiterations for each CFD time step. The Rayleigh time step is calculated based on the smallest particle size. DEM time steps are substantially

5. RESULTS AND DISCUSSION In this last section, the different parameters affecting the process will be successively discussed. After briefly comparing the results obtained by one-way and four-way coupling, the dynamics of isolated crystals will be considered through CFDDEM, before increasing the number of crystals in the reactor up to 10000. For all conditions, the impact of liquid volume and the occurrence of collisions will be quantified. In a final step, best seeding conditions are found by investigating uniform or localized seeding through a funnel for different crystal sizes. 5.1. Comparison of Predictions for One-Way and Four-Way Coupling. After the flow field has attained periodic steady state (5 s of physical time, corresponding to 25 revolutions of the impeller1), a single potassium dihydrogen phosphate (KDP) crystal with 0.5 mm size has been injected into the liquid. The location “R” defined in our previous investigation (see Figure 3 of ref 1) is considered as the injection point. This injection point is located close to the central axis in a region strongly affected by the impeller motion, 148



Article

Figure 3. Comparison of one-way and four-way coupling for (a) vertical position and (b) vertical velocity.

for T always corresponds to the total duration of the CFDDEM simulation (in physical time) and is up to 10 s. Two different injection locations are considered, point I of coordinates defined as (0.145,0.005,0.2 m) and point II of coordinates (0.05,0.005,0.2 m), shown in Figure 4a,c,

and it was observed that this location leads to unwanted sedimentation of crystals. It is therefore noteworthy to check whether this prediction could be modified when considering four-way coupling. Indeed, the obtained trajectories (not shown in the interest of saving space) are almost identical to that shown in ref 1, and sedimentation is observed as well in the CFD-DEM simulation, since the crystal density is noticeably higher than water density. Remembering that drag and buoyancy forces acting on particles are considered using similar models in both approaches, the obtained agreement is not unexpected for a single particle, as long as interactions are not important. In order to check quantitatively the differences between DPM and DEM, the temporal variation of the vertical position along the crystal trajectory and the magnitude of vertical velocity are compared in Figure 3a,b. A close agreement is observed. As expected, the sole difference between the oneway and the four-way simulation is the occurrence of wall collisions and the associated modeling. Here, collisions occur at the bottom of the reactor, appearing as visible oscillations in Figure 3. Even if they do not change the final outcome in the present case (sedimentation), they might lead to noticeable differences for other conditions, in which a crystal bouncing back from the bottom wall could be entrained again by the flow. 5.2. Single Crystal Dynamics. As aforementioned, the essential quantities needed to quantify local crystal growth in a crystallizer are the local slip velocity (velocity difference between particulate and liquid phase), the local crystal size distribution, and the local suspension density.30 In addition to its effect on crystal growth rate, the slip velocity will also affect classification of crystals. Therefore, we concentrate our investigation on this essential quantity in what follows. As a quantitative indicator proportional to crystal growth, the average slip velocity (ASV, in m/s) between liquid and crystal phase has been calculated based on the magnitude of the velocity difference as ASV =

1 no. of injected crystals

no. of injected crystal

∑ i=1

1 T

∫0

Figure 4. Crystal trajectories of single KDP crystal obtained for injection point I (top) and injection point II (bottom), for a liquid volume of 35 (left) or 24 L (right), colored by time.

respectively. Point I is chosen in the outer zone, outside of the circular baffle (draft tube), that is, in the zone that is less affected by the periodic motion of the impeller. Point II is placed closer to the central axis, that is, in a region strongly affected by the impeller motion. Remembering the importance of the liquid volume in the reactor observed in our previous study,1 the trajectories of single KDP crystals injected at points I and II are simulated through CFD-DEM for liquid volumes of 35 and 24 L, as shown in Figure 4 (left and right column, respectively). It is observed that the crystals are temporarily trapped by the

T

|Δvi(t )| dt

(14)

This is a quantity averaged in time over the total duration T. It will be used at first for a single crystal and then for a collection of up to 10000 crystals. As discussed above, this ASV influences the part related to mass transfer in the overall growth rate of crystals (see eq 13). In what follows, the value employed 149



Article

and (2) when the effect of the initial boundary condition has been sufficiently damped. In order to check this point, Figure 5 shows the instantaneous slip velocity in the 24 L crystallizer for both injection points during 10 s. Remember that the CFD-DEM simulation is started only after reaching the periodic steady state for the hydrodynamics through CFD, which takes 5 s. As a whole, the flow has therefore been simulated during 15 s or 75 impeller rotations. In conclusion, it is evident from Figure 5 that it is in principle sufficient to pursue the CFD-DEM simulation during 5 s of physical time in order to reach a periodic steady-state solution for the particles. For the present conditions, simulating 3 s of physical time would already be enough. Thus, even if this time is very short compared with real process times, it is sufficient to compute only around 25 impeller revolutions for hydrodynamics and additionally 25 impeller revolutions after adding the crystal particles in order to assess in a quantitative manner crystal dynamics in the reactor. It should be noted that computing 5 s of physical time might still take weeks or even months of computational time, depending on the employed computer and number of particles. It is, however, essential to involve many crystals in order to derive significant statistics and to investigate collisions in a realistic manner. The paramount importance of collisions in CFD-DEM simulations has already been discussed in scientific literature for different processes; see, for instance, ref 32. 5.3. Crystal Dynamics for up to 10000 Crystals. Now, the number of crystals injected in the CFD-DEM simulation is progressively increased from 1 up to 10000, considering again liquid volumes of 35 and 24 L. The instantaneous position of the crystals is shown in Figure 6 at the end of the CFD-DEM simulation. In the present case, a homogeneous injection has been implemented, implying that the initial crystal position is initiated randomly within the liquid volume when beginning the coupled simulation, corresponding initially to a perfectly mixed state. Confirming the previous results, it is apparent that many crystals are temporarily trapped by the rotating impeller when a crystallizer volume of 35 L is used (top row of Figure 6); however this phenomenon does not appear at a reduced liquid volume (24 L, bottom row). Since crystal-impeller and crystalbaffle collisions are indirectly responsible for secondary nucleation,27 the liquid volume of 24 L once again appears to be superior for practical purposes. It is also clear from the top row of Figure 6 that the upper part of the crystallizer does not participate in the process for a liquid volume of 35 L, since no crystals are found there at all. To investigate collision effects, the numbers of collisions of each crystal with the rotating impeller and with the stationary baffle counted during the entire simulation are reported in Table 2. In order to compare directly the results for widely varying crystal numbers, the total number of crystals in the reactor has been used to normalize the number of collisions. Table 2 contains engrossing information and requires extensive discussion. First, increasing the number of crystals from 1 to 10000 does not lead to a purely monotonic trend. Although such a trend still emerges, it is overlaid with random variations between the different cases. This is not necessarily surprising, since a random number generator is used initially to place the crystals in the liquid, thus impacting the collisions during the first time steps. Additionally, it might be necessary to take into account even more crystals before obtaining fully converged statistics.

rotational motion of the impeller for the liquid volume of 35 L (particularly visible in Figure 4a), an effect that does not appear for a crystallizer volume of 24 L. This first confirmation of our previous study is attributed again to the increase in average crystal velocity when the liquid volume decreases. In this sense, the lower liquid volume (24 L) is more favorable.1 However, sedimentation is an issue, since all four crystals in Figure 4 sediment within less than 10 s. Finding a suitable seeding position appears to be an essential issue. In order to now analyze crystal growth, the calculated values of slip velocity (ASV) are reported in Table 1. Since crystal growth is expected to increase linearly with ASV, as discussed previously, this quantity should be maximized.31 Table 1. Average Slip Velocity (in m/s, see eq 14) for Injection Points I and II, for Liquid Volumes of 35 and 24 L injection point I injection point II

liquid volume 35 L

liquid volume 24 L

0.0251 0.0203

0.0299 0.0205

It is clear from Table 1 that ASV is larger for the lower liquid volume (24 L). Additionally, the ASV strongly depends on the injection location (Figure 5), confirming previous findings

Figure 5. Evolution of the magnitude of slip velocity of single KDP crystal with time for injection points I and II with a liquid volume of 24 L.

concerning sedimentation.1 Crystals injected at point I (near the stationary domain) have a higher ASV than for those at point II. Finally, there is a 50% relative increase in ASV (and thus in crystal growth) when an injection at point II with 35 L of liquid is compared with an injection at point I in a 24 L volume. Though it is obviously too early to generalize (considering only isolated crystals and two locations), this first analysis reveals the very large impact of liquid volume and, even more, injection location on crystal growth. One important numerical issue is the needed duration of the simulation in physical time (parameter T). Incontestably, longer simulations lead to a higher computational effort, so shorter simulations would be advantageous. Conversely, a quantitative analysis based on a time average is only possible (1) after reaching a periodic steady state for flow and particles 150



Article

Figure 6. Crystal dynamics in crystallizer for a liquid volume of 35 L (a−e, top row) or 24 L (h−l, bottom row), considering an increasing number of crystals, from 100 to 10000 from left to right (100, a, h; 500, b, i; 1000, c, j; 2500, d, k; 10000, e, l). The results shown correspond to the end of the simulation, and crystals are colored by instantaneous velocity magnitude.

collisions with the baffle dominate (roughly 4 times more often than with the impeller). In order to quantify crystal growth, the average slip velocity between the crystal and fluid is calculated again using eq 14. Table 3 depicts the calculated average slip velocity (ASV), as

Table 2. Number of Collisions Per Crystal with the Rotating Impeller and with the Stationary Baffle for Liquid Volumes of 35 and 24 L seeding crystal numbers

at impeller (35 L)

at baffle (35 L)

at impeller (24 L)

at baffle (24 L)

1 100 500 1000 2500 10000

40 332.42 379.79 405.39 412.93 386.52

3 85.79 54.11 29.35 117.27 147.74

9 1.81 2.01 17.22 7.42 6.57

0 16.33 12.79 15.69 17.06 28.87

Table 3. Average Slip Velocity (ASV, see eq 14) in Liquid as Well as Average Relative Velocity (RV) When Crystals Collide with the Impeller for Liquid Volumes of 35 and 24 L

Bearing in mind that the simulation with 10000 crystals leads to approximately one month of computing time, testing this issue will be part of our future work. Since the statistical content of the simulation increases with the number of crystals, all values discussed next correspond to the case with 10000 crystals (last line in Table 2). The scatter observed in the data for 35 L liquid volume is probably due to the large number of particles trapped within the impeller region and sedimenting at the bottom of the reactor. For this reason, the statistics for the 24 L volume are more compelling and also more relevant for practical purposes. However, the qualitative trends are the same for both liquid volumes. It is particularly evident from Table 2 that the average number of collisions per crystal hitting the impeller or stationary baffle is much higher for the larger liquid volume of 35 L. The crystal movement is slower in this case, and many crystals are trapped by the impeller, as previously discussed. For a liquid volume of 24 L, there are roughly 60 times fewer collisions with the impeller and 5 times fewer collisions with the baffle. This information is crucial when investigating secondary nucleation rate in crystallizers, because it is proportional to the product of collision energy and frequency of collision.33 While for the large liquid volume, and thus for the slower crystals, collisions occur mostly with the impeller (roughly 3 times more than with the baffle), the situation reverses for the lower liquid volume and faster crystals, where

seeding crystal numbers

ASV (35 L) (m/s)

ASV (24 L) (m/s)

RV (35 L) (m/s)

RV (24 L) (m/s)

1 100 500 1000 2500 10000

0.0203 0.0238 0.0233 0.0241 0.0255 0.0286

0.0205 0.0326 0.0329 0.0338 0.0358 0.0401

1.20 1.50 16.81 20.02 21.04 25.23

0.26 0.42 0.55 1.04 1.16 1.24

well as the relative velocities (RV) observed when crystals collide with the impeller. All data confirms the advantage of the reduced liquid volume, 24 L. First, the slip velocity (ASV) is systematically higher for the 24 L reactor, 40% higher when averaging over 10000 crystals, indicating a similar difference in crystal growth rate. At the same time, the relative velocity during collisions is much smaller with a liquid volume of 24 L, typically by a factor 20. This is because, as aforementioned, crystal dynamics and particle velocity distribution are completely different in the 24 L and in the 35 L reactor; the crystals moving in the latter do so at a much slower speed. Consequently, crystal attrition and secondary nucleation are expected to be much higher for a liquid volume of 35 L. When a small number of particles is considered, a large proportion of those appear to be trapped within the plane of the impeller, in particular for a liquid volume of 35 L. When 500 crystals or more are considered, this effect becomes progressively less important. This is partly due to increased interactions ejecting the crystals outside of the impeller plane. Now, it might be useful to compute the total energy loss, which is the sum of normal and tangential energy loss during 151



Article

5.5. Uniform vs Localized Seeding. Seeding has become one of the most critical steps for optimizing process efficiency and product quality and for ensuring the final particle size. Inconsistent PSD can often be traced back to inappropriate seeding. Many parameters must be taken into consideration when designing a seeding strategy such as seed size, seed loading, and seed addition temperature. These parameters are generally optimized manually based on process kinetics and the desired final particle properties. In a real crystallization process, the seeding of crystals is usually carried out using funnel injection. To understand such seeding effects, 1000 crystals are now injected in the CFD-DEM simulation through a realistic funnel, similar to that used in real operation, modeled with an upper radius of 2 cm and an end radius of 1 cm, with a height of 1.7 cm. Three different seeding locations are compared as shown in Figure 7: funnel position I (SD, near the reactor wall) corresponds to (x, y, z) = (0.05, 0.005, 0.25 m); funnel position II (middle, above the baffle) corresponds to (0.092, 0.005, 0.25 m); funnel position III (RD, near the rotation axis) is placed at (0.133, 0.005, 0.25 m). These coordinates correspond to the center of the lower boundary of the funnel. In all three cases, crystals do not appear to be trapped by the rotation of the impeller, confirming previous statements for the lower liquid volume. To quantitatively optimize the funnel location, two conditions are once again crucial: maximizing ASV and minimizing collisions. Both parameters are listed in Table 6 for the three different funnel positions. A clear optimum appears from Table 6. It is found that the number of collisions per crystal is minimum when the funnel is in position II, above the stationary baffle. At the same time, this position leads to the highest ASV value, almost three times better than the poorest choice (injecting near the rotation axis, position III). Therefore, seeding using a funnel above the stationary baffle is recommended. Comparing now Table 6 to Table 3, the ASV value is doubled when using a funnel in optimal position instead of injecting the seeds homogeneously within the entire reactor, therefore confirming the advantage of the technical solution used in practice. 5.6. Effect of Crystal Size. Finally, the impact of seed crystal size on the findings is investigated. For this purpose, crystals of different sizes are injected simultaneously in different amounts from funnel position II (above stationary baffle), previously found to be the best solution. Figure 8a shows the final distribution of the crystals within the crystallizer. The corresponding crystal size distribution initiated in the funnel is shown in Figure 8b and corresponds to a normal (Gaussian) distribution in terms of crystal volume (or mass). The resulting ASV and the number of collisions per crystal at the rotating impeller are reported in Table 7 as a function of the size of the seed crystals. Table 7 indicates that there is a crystal size for which a minimum number of collisions per crystal are obtained. For a crystal diameter from 0.1 to 0.5 mm, the average number of collisions per crystal decreases rapidly, before increasing again slowly for larger crystals. This indicates that crystals around 0.5 mm in dimension would be less subjected to attrition. The average slip velocity (ASV) increases slightly with crystal diameter, but the variation is quite small, about 5% relative change for crystals between 0.1 and 1 mm in diameter. Hence, the ASV can be assumed to be quite independent from crystal

the collision process for 35 and 24 L, see Table 4. The previous results confirm that the collision energy loss per crystal for the Table 4. Total Energy Loss per Crystal Due to Collisions of Crystal with the Rotating Impeller and with the Stationary Baffle for Liquid Volumes of 35 and 24 L seeding crystal numbers 1 100 500 1000 2500 10000

total energy loss per crystal (J) (35 L) 0 2.98 1.93 2.33 1.48 7.79

× × × × ×

10−12 10−12 10−7 10−7 10−8

Total energy loss per crystal (J) (24 L) 0 8.58 1.96 1.03 4.56 1.77

× × × × ×

10−13 10−13 10−13 10−14 10−14

smaller liquid volume (24 L) is orders of magnitude smaller than for the larger volume (35 L). For the higher liquid volume, energy losses are dominated by intense collisions with the impeller, leading to high secondary nucleation rates30 for these conditions. In conclusion, the crystal growth is found to be faster when a liquid volume of 24 L is used. The number of collisions, relative velocity of those, and total energy loss associated with collisions are considerably lower for 24 L. Therefore, the CFD-DEM simulation indicates that a 24 L liquid volume is highly advantageous in every respect. Consequently, further investigations focus on the 24 L crystallizer volume. 5.4. Impeller Rotation Rate. The rotational speed of the impeller plays an indisputably significant role as a process parameter and may impact the collision number, collision intensity, and crystal growth in a complex manner. Concentrating now on the liquid volume of 24 L, the simulations have been repeated for three different impeller speeds: 250, 300, and 350 rpm. Only 1000 crystals are involved in the CFD-DEM simulations in order to speed-up the computations. Resulting data for both average slip velocity and number of collisions per crystal can be found in Table 5. When we analyze this data, the Table 5. Average Slip Velocity and Number of Collisions Per Crystal with the Rotating Impeller As a Function of Impeller Rotational Speed for 24 L Liquid Volume rotational speed of impeller (rpm)

ASV (m/s)

collisions/crystal at impeller

250 300 350

0.0232 0.0338 0.0448

5.36 17.22 29.02

expected outcome that both ASV and average number of collisions increase monotonically with impeller rotational speed is observed. In particular, the ASV appears to increase proportionally to the square of the rotational speed of the impeller. Therefore, the ASV is proportional to the energy input. However, this correlation is only based on three values and must be checked more extensively in the future. The number of collisions with the impeller per crystal increases even faster with the rotational speed. Therefore, considering that an operator tries to maximize ASV (hence, crystal growth) while minimizing secondary nucleation and crystal attrition, an optimal impeller speed will be found in the intermediate range, confirming practical experience: neither very low nor very high impeller speeds are appropriate. 152



Article

Figure 7. Crystal dynamics in crystallizer for a liquid volume of 24 L considering 1000 crystals injected through a funnel placed at three different positions. Results correspond to the end of the simulation and crystals are colored by instantaneous velocity magnitude.

Table 6. Average Slip Velocity and Number of Collisions Per Crystal with the Rotating Impeller As a Function of the Funnel Position Used for Injection funnel location

ASV (m/s)

collision/crystal (impeller)

position I position II position III

0.0397 0.0716 0.0254

3.221 1.189 1.834

Table 7. Average Slip Velocity and Number of Collisions per Crystal with the Rotating Impeller as a Function of Crystal Diameter for a Gaussian CSD at Injection

size for the considered conditions, so that the evolution of the CSD would be mostly constrained by attrition. These observations do not change much when we switch from a Gaussian CSD to a top-hat CSD (compare Tables 7 and 8, in which the same number of crystals has been injected through the funnel for each size class). Here again, crystals with a diameter of 0.5 mm are associated with a minimum number of collisions. The relative change in ASV is below 9% for a factor 10 in crystal diameter. Therefore, it is acceptable to consider the ASV as constant over this range. The observed differences in collision probability are far more considerable. As a trade-off between crystal growth and crystal attrition, KDP

seeding crystal number

size of the crystal (mm)

ASV (m/s)

collision/crystal (at impeller)

50 200 500 200 50

0.1 0.3 0.5 0.7 1

0.0278 0.0281 0.0284 0.0301 0.0305

2.04 1.38 0.42 0.62 0.64

seed crystals with a diameter of 0.5 mm can therefore be considered as optimal for the designated process.

6. CONCLUSIONS The one-way coupling (DPM) is based on a translational force balance that is formulated for individual particles or parcels of particles. As a major difference from DEM, particle collisions are usually neglected with DPM and the rotation of particles is not considered. Additionally, the volume fraction of the particulate phase is not considered in DPM. Due to these assumptions and simplifications, DPM is only valid for dilute

Figure 8. (a) Crystal dynamics in crystallizer for a liquid volume of 24 L considering 1000 crystals of different sizes injected through the funnel in position II. Results correspond to the end of the simulation and crystals are colored by instantaneous velocity magnitude. (b) Initial crystal size distribution within the funnel. 153



Article

Present Addresses

Table 8. Average Slip Velocity (in m/s) and Number of Collisions Per Crystal with the Rotating Impeller as a Function of Crystal Diameter for a Top-Hat CSD at Injection seeding crystal number

size of the crystal (mm)

ASV (m/s)

collision/crystal (at impeller)

200 200 200 200 200

0.1 0.3 0.5 0.7 1

0.0281 0.0282 0.0284 0.0297 0.0295

1.28 1.13 0.86 1.04 1.27

∥

B.A.A.: National Institute of Technology Karnataka, Mangalore, India. ⊥ M.B.: Hüttlin GmbH, a Bosch Packaging Technology Company, Schopfheim, Germany. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Dr. B. Ashraf Ali would like to sincerely thank Pro3 (http:// www.verfahrenstechnik-pro3.de), the German Process Engineering Expertise Network, for the financial support of his work. Additionally, the support of the Deutsche Forschungsgemeinschaft (DFG) within the Research Program SPP 1679 “Dynamische Simulation vernetzter Feststoffprozesse” is gratefully acknowledged. Fruitful discussions with Mrs. Kristin Kerst and Mr. Erik Temmel have helped improve the paper. Concerning EDEM, the support from Mr. Senthil Arumugam (Academic Program Manager), DEM Solutions Ltd, as well as his team are gratefully acknowledged.

fluid−particle flows for which collisions are not essential. The main advantage over DEM is, however, that large time steps can be used, leading to much smaller computational times. As an alternative, fully coupled CFD-DEM simulations, as considered in this work, require a significant amount of computational effort in terms of computer performance and memory. This is the only way to take into account in a realistic manner collision processes and locally dense conditions. The crystal dynamics in an existing batch crystallizer have been investigated using four-way coupling by injecting first single KH2PO4 (KDP) crystals into water in CFD-DEM simulations. Though trajectories of isolated crystals are nearly the same, DEM allows accounting for collisions with rotating impeller, baffles, and reactor walls. Considering that collisions are essential to understand attrition and secondary nucleation, the CFD-DEM methodology should be preferred when possible. Up to 10000 individual crystals have been considered in the simulation. The impact of various process parameters on crystal growth has been quantified by analyzing the average slip velocity and collision rates. One notable discovery is that a liquid volume of 24 L is solely advantageous compared with a larger liquid volume of 35 L. When the crystals are injected in a realistic manner through a funnel, the obtained results strongly depend on the funnel position. The best injection location is found above the stationary baffle. While increasing impeller rotational speed is favorable for crystal growth, it simultaneously leads to rapidly increasing collision rates and thus attrition. A crystal diameter of 0.5 mm is associated with the minimum number of collisions and would therefore contribute mostly to crystal growth for the present process. To summarize the recommendations resulting from this study: a relatively low liquid volume should be used; injection should be accomplished from a funnel above the stationary baffle; injecting crystals of 0.5 mm leads to less collisions; and the rotation rate should be selected at an intermediate level, supporting the value of 300 rpm found experimentally.

■

ASSOCIATED CONTENT

■

AUTHOR INFORMATION

■

REFERENCES

(1) Ashraf Ali, B.; Janiga, G.; Temmel, E.; Seidel-Morgenstern, A.; Thévenin, D. Numerical analysis of hydrodynamics and crystal motion in a batch crystallizer. J. Cryst. Growth 2013, 372, 219−229. (2) Mullin, J. W.: Crystallization; Butterworth-Heinemann, Oxford, 2001. (3) Rohani, S.; Horne, S.; Murthy, K. Control of product quality in batch crystallization of pharmaceuticals and fine chemicals. Part 1: Design of the crystallization process and the effect of solvent. Org. Process Res. Dev. 2005, 9 (6), 858−872. (4) Rane, C. V.; Ganguli, A. A.; Kalekudithi, E.; Patil, R. N.; Joshi, J. B.; Ramkrishna, D. CFD simulation and comparison of industrial crystallizers. Can. J. Chem. Eng. 2014, 9999, 1−19. (5) Myerson, A. S.: Handbook of Industrial Crystallization; Butterworth-Heinemann, 2002. (6) Randolph, A. D. and Larson, M. A.: Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization, 2nd ed., Academic Press, New York, 1988. (7) Moscosa-Santillán, M.; Bals, O.; Fauduet, H.; Porte, C.; Delacroix, A. Study of batch crystallization and determination of an alternative temperature-time profile by on-line turbidity analysis application to glycine crystallization. Chem. Eng. Sci. 2000, 55 (18), 3759−3770. (8) Garside, J.; Gibilaro, L. G.; Tavare, N. S. Evaluation of crystal growth kinetics from a desupersaturation curve using initial derivatives. Chem. Eng. Sci. 1982, 37 (11), 1625−1628. (9) Doki, N.; Kubota, N.; Sato, A.; Yokota, M.; Hamada, O.; Masumi, F. Scaleup experiments on seeded batch cooling crystallization of potassium alum. AIChE J. 1999, 45 (12), 2527−2533. (10) Bohlin, M.; Rasmuson, A. C. Application of controlled cooling and seeding in batch crystallization. Can. J. Chem. Eng. 1992, 70 (1), 120−126. (11) Doki, N.; Kubota, N.; Sato, A.; Yokota, M. Effect of cooling mode on product crystal size in seeded batch crystallization of potassium alum. Chem. Eng. J. 2001, 81 (1−3), 313−316. (12) Lewiner, F.; Févotte, G.; Klein, J. P.; Puel, F. Improving batch cooling seeded crystallization of an organic weed-killer using on-line ATR FTIR measurement of supersaturation. J. Cryst. Growth 2001, 226 (2−3), 348−362. (13) Bałdyga, J. and Bourne, J. R.: Turbulent Mixing and Chemical Reactions; Wiley, 1999. (14) Marchisio, D. and Fox, R. O.: Computational Models for Polydisperse Particulate and Multiphase Systems; Cambridge University Press, 2013.

S Supporting Information *

A video file showing the time-dependent evolution of 10,000 crystals in the 24 L reactor; for this simulation, the crystals are initially homogeneously distributed in the crystallizer in a random manner. This material is available free of charge via the Internet at http://pubs.acs.org. Corresponding Author

*E-mail: [email protected]. 154



Article

(15) Bordás, R.; John, V.; Schmeyer, E.; Thévenin, D. Numerical methods for the simulation of an aggregation-driven droplet size distribution. Theor. Comp. Fluid Dyn. 2013, 27, 253−271. (16) Wei, H.; Zhou, W.; Garside, J. Computational fluid dynamics modeling of the precipitation process in a semibatch crystallizer. Ind. Eng. Chem. Res. 2001, 40 (23), 5255−5261. (17) Liiri, M.; Hatakka, H.; Kallas, J.; Aittamaa, J.; Alopaeus, V. Modelling of crystal growth of KDP in a 100 dm3 suspension crystallizer using combination of CFD and multiblock model. Chem. Eng. Res. Des. 2010, 88 (9), 1297−1303. (18) Rane, C. V.; Ekambara, K.; Joshi, J. B.; Ramkrishna, D. Effect of impeller design and power consumption on crystal size distribution. AIChE J. 2014, 60 (10), 3596−3613. (19) Zauner, R.; Jones, A. G. Scale-up of continuous and semibatch precipitation processes. Ind. Eng. Chem. Res. 2000, 39 (7), 2392−2403. (20) Jaworski, Z.; Nienow, A. W. A strategy for solving a mathematical model of barium sulphate precipitation in a continuous stirred tank reactor. Chem. Eng. Res. Des. 2004, 82 (9), 1089−1094. (21) Tsuji, Y.; Tanaka, T.; Ishida, T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol. 1992, 71 (3), 239−250. (22) Hobbs, A. Simulation of an aggregate dryer using coupled CFD and DEM methods. Int. J. Comput. Fluid Dyn. 2009, 23 (2), 199−207. (23) Deen, N. G.; Van Sint Annaland, M.; Van der Hoef, M. A.; Kuipers, J. A. M. Review of discrete particle modeling of fluidized beds. Chem. Eng. Sci. 2007, 62, 28−44. (24) EDEM2.1: Theory and Reference guide; DEM solutions Ltd, Edinburgh, 2010. (25) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci. 2007, 62, 3378−3396. (26) Mersmann, A., Eble, A. and Heyer, C.: Crystallization Technology Handbook; Marcel Dekker Inc., New York, 2001. (27) Liiri, M.; Koiranen, T.; Aittamaa, J. Secondary nucleation due to crystal-impeller and crystal-vessel collisions by population balances in CFD-modelling. J. Cryst. Growth 2002, 237−239 (1−4 III), 2188− 2193. (28) Mersmann, A.; Braun, B.; Löffelmann, M. Prediction of crystallization coefficients of the population balance. Chem. Eng. Sci. 2002, 57 (20), 4267−4275. (29) Kim, S.; Myerson, A. S. Metastable solution thermodynamic properties and crystal growth kinetics. Ind. Eng. Chem. Res. 1996, 35 (4), 1078−1084. (30) Hatakka, H., Liiri, M., Aittamaa, J., Alopaeus, V., Kultanen, L. M. and Kallas, J.: Flow patterns and slip velocities of crystals in a 100-liter suspension crystallizer equipped with two turbine impellers, In Proc. 17th International Symposium on Industrial Crystallization - ISIC 17 CGOM8, 2008, 1851−1858. (31) Liiri, M.; Enqvist, Y.; Kallas, J.; Aittamaa, J. CFD modelling of single crystal growth of potassium dihydrogen phosphate (KDP) from binary water solution at 30 °C. J. Cryst. Growth 2006, 286 (2), 413− 423. (32) Fries, L.; Antonyuk, S.; Heinrich, S.; Palzer, S. DEM-CFD modeling of a fluidized bed spray granulator. Chem. Eng. Sci. 2011, 66 (11), 2340−2355. (33) Ottens, E. P. K.; Janse, A. H.; De Jong, E. J. Secondary nucleation in a stirred vessel cooling crystallizer. J. Cryst. Growth 1972, 13−14 (C), 500−505.

155


Coupled Computational Fluid DynamicsâDiscrete Element Method

Recommend Documents

Coupled Computational Fluid DynamicsâDiscrete Element Method