Mathematical Modeling of Layer Crystallization on a Cold Column

Article pubs.acs.org/IECR

Mathematical Modeling of Layer Crystallization on a Cold Column with Recirculation Nima Yazdanpanah, Allan Myerson, and Bernhardt Trout* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States Novartis-MIT Center for Continuous Manufacturing, Cambridge, Massachusetts 02139, United States S Supporting Information *

ABSTRACT: Falling film crystallization is a novel technique used to overcome the challenges of solid and slurry handling used for continuous manufacturing in the pharmaceutical industry. The process and its performance were demonstrated recently. In this study, a dynamic model of falling film crystallization with a recirculation loop was developed in the form of coupled partial differential equations describing the effects of mass transfer, heat transfer, impurity inclusion, and crystallization kinetics. The model described the variation of processing parameters, such as the flow rate and cooling temperature, for three different crystallization solution systems (16 cases). The experimental results confirm the predictions of the model, which will provide further insights into processes in which experimental values do not exist or are difficult to obtain, such as the axial temperature profile of the falling liquid film or the overall purity of the deposited crystal layer at different time intervals.

1. INTRODUCTION Crystallization is a common purification and separation technique used in the pharmaceutical industry to separate the drug product from the liquid phase, which also contains impurities. However, handling of the slurry solution, filtration, and washing and drying the crystals are challenging for batch and continuous manufacturing. Fouling and clogging of the lines in both continuous crystallization and flow chemistry methods still present a challenge for designing continuous processes. The filtration, washing, and drying of the crystals reduce the overall yield of the process, and poor filtering causes many issues in downstream processes.1 These challenges have prompted the development of a novel technique to avoid solid handling and filtration for both batch and continuous manufacturing. The other objective is to maintain, and even improve, the overall yield and purity of the final product, which are the key efficiency parameters. The novel technique of falling film solution crystallization has recently been developed to eliminate the need for solid handling and filtration operations by directly crystallizing the main compound in solution on a coldfinger.2 The deposited crystal layer on the coldfinger, which has a high level of purity, then dissolves in a fresh solvent and is transferred for downstream processing and formulation. The process parameters, such as the flow rate, temperature, and solution system (solubility profile), play significant roles in the final purity of the product, overall yield of the system, and the performance and control of the system. Although experimental results published previously2 demonstrate the overall system perform© 2016 American Chemical Society

ance for a range of crystallization solution systems and process parameters, it is necessary to develop a mathematical model for further improvement of the falling film crystallizer and to predict the performance for new crystallization cases. The falling film crystallization technique is highly effective due to the rejection of impurities on the growing crystal layer interface, which are transported away from the crystallizer. Therefore, the growth rate (temperature and concentration dependent) and mass transfer (through the flow rate) significantly affect the process. The solution used in the falling film crystallizer is very similar to that used in melt crystallization, although the influencing parameters are different; in melt crystallization, the role of the heat transfer rate is more important, where the mass transfer rate is more governing for solution crystallization.3−5 In the recent years, many modeling studies of melt crystallization have been performed.4,6−13 However, because falling film solution crystallization is a new technique, no mathematical models have been developed for this method. For layer melt crystallization, the experimental and modeling research on crystal layer growth has focused on how the operational conditions affect the process efficiency and crystal growth. Some attempts have been reported to develop models of the impurity distribution, mostly based on theoretical calculations using a solid−liquid phase Received: Revised: Accepted: Published: 5019

January 13, 2016 March 27, 2016 April 11, 2016 April 20, 2016 DOI: 10.1021/acs.iecr.6b00057 Ind. Eng. Chem. Res. 2016, 55, 5019−5029

Article

Industrial & Engineering Chemistry Research

Figure 1. Schematic description of state variables for the column.

diagram of the mixture, and of the dependency of practical purification on many critical operational conditions, such as the crystal growth rate and transport phenomena.12−21 In general, the few reported models have attempted to predict the overall yield and purity of the product of layer melt crystallization. The aim of this work is to develop a mathematical model for the falling film solution crystallizer and to predict the performance of this process, in recirculation mode, under various operation conditions for three different crystallization solution systems. The numerical results are compared to and validated by the experimental results to predict the crystal thickness (yield) and final purity of the product.

6. There is no heat transfer to the surrounding air. The heat transfer occurs by conduction with the coolant in the column through the pipe wall and deposited crystal layer. Heat transfer for the liquid side (falling liquid film) occurs by convection. 7. There is no evaporation from the surface of the film. 8. The feed tank is well-mixed (uniform temperature and concentration) in the recirculation mode. 9. The nucleation effect is negligible. Crystal growth occurs on the surface of the column and deposited crystals, and no nucleation and growth occurs in the bulk solution. 10. The latent heat of crystallization is negligible in comparison with the heat flux from cooling of the falling solution film. 2.2. Gridding and Assignment of State Variables for Differential Elements. The schematic diagram in Figure 1 shows the assignment of state variables for the internally cooled column surrounded by a deposited layer of crystals and a falling liquid film. The direction of heat flux is from the center of the column toward the surfaces, and the direction of mass transfer is along the reverse path. The rm and Tm are constants at the surface of the column (metal pipe) and the rs,Ts and rf,Tf change over time due to the growing layer of crystals. The length of the column, L, is divided equally into a number of difference elements, i, in the axial direction, Δxi. 2.3. Mass Transfer. 2.3.1. Mass Transfer of Solute in the Liquid Phase. According to the above assumptions, the solute mass balance equation is represented by Fick’s diffusion law in cylindrical coordinates as follows:

2. MODEL DESCRIPTION A model of the continuous crystallization process based on the simultaneous solution of the heat and mass balance equations was developed. The core of the model is similar to our recently published work;22 however, the models for growth rate and kinetics estimation and calculation of impurity inclusion and purification are different. The growth rate section of the model is based on the mass deposition rate over the column and the growth rate parameters calculated from the desupersaturation profiles, accordingly. The distribution coefficients for different temperatures, impurity contents, and process kinetic parameters are incorporated into the model to obtain the final crystal purity. The purpose of the model is to predict the crystal purity, yield, and best processing time for the continuous crystallization system on a cold column, and the results will be used for further development of the process using different numbers of stages and processing parameters. Hence, the fitted crystal growth parameters are the best guess values; however, the model is still a useful tool to investigate the purity and yield of the product. 2.1. Assumptions. The following are several assumptions and limitations of the model: 1. The cooling temperature of the coolant through the column is constant, due to the high mass flow rate of cooling water and adequate cooling duty available to cause negligible temperature rise for the cooling water. 2. The solution flow over the column consists of a fully developed laminar flow. The interface wave effect and interfacial shear are negligible. 3. No concentration profile is present in the thin film of the solution in the radial direction. The solution in the differential volume of the film (over a length scale) is a homogeneous mix of the solvent, solute, and impurities. 4. The concentration profile of impurities at the interface of the crystals is negligible, and the solution is well-mixed. 5. The thermo-physical properties of the solution are independent of temperature.

1 ∂ ⎛⎜ ∂C ⎞⎟ ∂C ∂C r , + u(r ) =D r ∂r ⎝ ∂r ⎠ ∂t ∂x

rs ≤ r ≤ r f ,

0≤x≤L (1)

where C is the solute concentration, u(r) is the fluid velocity, D is the molecular diffusion coefficient, x is the axial distance, r is the radial distance, and L is the total length of the column. u(r), the velocity for fully developed laminar follow for a falling film in cylindrical coordination, is defined by eq 2. u(r ) =

ρg ⎛ 2 r⎞ ⎜rs − r 2 + 2r f 2 ln ⎟ , 4μ ⎝ rs ⎠

rs ≤ r ≤ r f

(2)

At the solid−liquid interface, the mass flux to the surface is equal to the crystallization rate. Therefore, the initial and boundary conditions for the mass transfer equation for the solute in liquid phase are 5020

DOI: 10.1021/acs.iecr.6b00057 Ind. Eng. Chem. Res. 2016, 55, 5019−5029

Article


phases at equilibrium. Therefore, in the impurity mass balance equation, the initial and boundary conditions are similar to those of the solute mass balance, which is defined by eq 6.

IC: t = 0, C(r , x , 0) = C0 x = 0,

∂Cimp , sol

dC(r , 0, t ) dt F(L , t )(Csol(L , t ) − Csol(0, t )) = Vtank

∂t

rs ≤ r ≤ r f ,

∂Cimp , sol ∂x

= Dimp , sol

1 ∂ ⎛ ∂Cimp , sol ⎞ ⎜r ⎟, r ∂r ⎝ ∂r ⎠

0≤x≤L

(6)

IC:

BC:

t = 0, Cimp , sol(r , x , 0) = Cimp , sol ,0

∂C r = r f , −D =0 ∂r

x = 0,

∂C r = rs , −D = −G = −ki(Cr − Csat )n ∂r

r = rs , − Dimp , sol

r

(3)

The solute mass flux, j, and mass transfer coefficient, kd, at the crystal−solution interface are defined by eqs 4 and 5. js , r = −D

∂C ∂r

js , r Cr − Csol

(4)

D ShD ∂C =− = Cr − Csol ∂r z

∂Cimp , sol ∂r

∂Cimp , sol ∂r

=0

⎛ ⎞ Cimp , sol = − G⎜⎜Cimp , sol − K ρH , s ⎟⎟ CH , sol ⎝ ⎠

Cimp,sol is the concentration of impurities in the solution, K is the distribution coefficient, CH,sol is the concentration of the host compound (main solute compound) in the solution, and ρH,s is the density of the host compound in the solid phase. Distribution Coefficients. Distribution coefficient is defined as the ratio of the impurity concentration to the host compound concentration in the solid phase divided by that ratio in the liquid phase. The following equation defines the distribution coefficient:

s

s

Vtank

r = r f , − Dimp , sol

∫r f C(r , x , t )U (r , t )r dr ∫r f u(r , t )r dr

dt F(L , t )(Cimp , sol(L , t ) − Cimp , sol(0, t ))

BC:

r

Csol(x , t ) =

dCimp , sol(r , 0, t ) =

where F is the flow rate of the solution, C0 is the initial concentration, Csol is the average solution concentration at distance x, G is crystal growth rate, ki is the rate of solute integration during crystal formation, Cr is the concentration at the crystal layer surface, Csat is the concentration at the saturation point (solubility), and n is the kinetic exponent rate for crystal growth. Csol, the average concentration, is defined by eq 3 and calculated by measuring the outlet concentration from the bottom of the column and the desupersaturation curve.

kd =

+ u(r )

Distribution Coefficient =

(Cimp/CH )s , final (Cimp/CH )sol ,0

(7)

where Cimp is the concentration of the impurities and CH is the concentration of the host compound. Two types of distribution coefficients mainly exist, “equilibrium” (Keq) and “effective” (K) coefficients. The equilibrium distribution coefficient of an impurity in a solute is a function of the solution thermodynamics and solid-state thermodynamics. The effective distribution coefficient is regulated by the process parameters, such as the cooling rate, mixing, local concentration of the impurity, solvent inclusion, and crystal growth rate. Therefore, the effective distribution coefficient is the ratio of the impurity concentration in the crystal to its average concentration at the crystal−medium interface. This concentration at the interface is different from that in the bulk medium because the growth interface acts either as a source or as a sink. The effective conditions that result in the inclusion of more impurities in the crystals in the form of surface impurities, lattice impurities, and solvent inclusion increase the distribution coefficient to an effective level. The contribution of the kinetic effects on the purity of the crystals or crystal layers will depend on the configuration of the crystallizer and the process conditions during crystal growth. The concentration of impurities at the interface may be significantly increased due to mass transfer, increasing the tendency to incorporate impurities. This produces a strong contribution to the overall distribution coefficient when a large growth rate is established by elevated supersaturation/supercooling. Ulrich and Neumann24 measured the effective

(5)

Figure 2 gives a schematic description of transport phenomena and changes in concentrations and fluid temperature in the axial position along the length of the column. 2.3.2. Mass Transfer of Impurities. The measurement of the inclusion of impurities in the crystal lattice has been described previously in the form of distribution coefficients.23 The measured distribution coefficient at different impurity levels, K, is used for the segregation of impurities in the liquid and solid

Figure 2. (A) Schematic description of heat transfer and profile of the temperature. (B) Changes in concentrations and fluid temperature in the axial position along the length of the column. 5021


Article

Industrial & Engineering Chemistry Research distribution coefficient in which the effect of flow (operating conditions) on the effective distribution coefficient was demonstrated. Jiang et al.3 proposed an additive type effective (average) distribution coefficient function that considers impurity trapping as being dependent upon the crystal network type/structure. Considering the kinetics and thermodynamic effects, an average distribution coefficient equation was used, which consisted of a combination of the equilibrium distribution coefficient, impurities from solvent inclusion, and entrapment of impurities due to crystal growth. ⎛ msolinc , t ⎞ × Cimp , t ⎟⎟ + KG K = (Keq(RelCimpt )) + ⎜⎜ ⎝ ms , t ⎠

therefore: G=

(8)

∂T ∂T ∂ ⎛ ∂T ⎞ + u(r ) = α ⎜r ⎟ , ∂t ∂x ∂r ⎝ ∂r ⎠

dt

t = 0, T (r , x , 0) = T0 x = 0, T (r , 0, t ) = Ttank

BC:

δ ΔCimp , tG exp −

= ki(Cr − Csat )n A =

GδD D

)

r = r f , −ksol

(9)

∫ Gρs dA

r = rs , −ksol

∂T =0 ∂r

∂T = Ucol(Tcol − Ts) ∂r

where T is the temperature, α is the thermal diffusivity, ksol is the thermal conductivity of the solution, Tcol is the cooling temperature of the column, and Ucol is the overall heat transfer coefficient from eq 14. ⎛ 1 ln(rm/rcol ,0) ln(rs/rm) ⎞ 1 ⎟⎟ = ⎜⎜ + + Ucolrs km ks ⎠ ⎝ hwrcol , i

(14)

km and ks are the thermal conductivities of the metal and solid, respectively, hw is the heat transfer coefficient of the cooling water in the column, and rcol is the radius of the column.

3. NUMERICAL METHOD The equations in the model were solved using the ODE23 solver in Matlab. The method-of-line technique was used to discretize the domain and solve the equation for each finite difference element. The interchange coefficients were determined, and the heat and mass transfer governing equations were solved according to the following steps: 1. Assuming the initial and boundary conditions (Tables S.2 and S.3). 2. Calculating the growth rate parameters over a time scale. 3. Solving the heat and mass transfer equations simultaneously for the first length step. 4. Updating the fluid temperature and concentration for the next length step. 5. Defining the heat and mass fluxes simultaneously for the total length.

(10)

where dmH,s/dt is the rate at which the solute mass is transported over the area, A = 2πrsΔx. Therefore, the total mass in the differential length element of Δx in the axial direction is Δṁ = ρπ (rs 2 − rm 2)Δx s

(13)

IC:

D is the diffusivity, ΔCH,T is the change in concentration of the host compound at a given time, δ is the thickness of the liquid film, ΔCimp,t is the change in concentration of impurities at a given time, G is the crystal growth rate, and δD is the thickness of the mass transfer boundary layer. The key parameters for determining the behavior of impurities in the growth of the crystalline layer are the temperature and concentration profiles of the boundary layer near the solid−liquid interface. It is difficult to predict the local values of the effective distribution coefficient and the time changes in the concentration and temperature profiles at the boundary layer and the solid−liquid interface. However, the interfacial distribution coefficient varies with the growth of the crystalline layer and the elapsed time. Therefore, the value of KG was regressed from the change in concentration of the impurities from the experimental results. 2.3.3. Crystal Growth (Desupersaturation). The desupersaturation profiles for different temperatures and flow rates were generated by measuring the change in concentration of the active pharmaceutical ingredient (API) in the feed tank in the experimental results. These profiles are the basis for measuring the growth rate, G, and the rate of solute integration in the formation of crystals, ki. The growth rate of the crystal layer is derived from the mass balance at the solid−liquid interface: dmH , s

rs ≤ r ≤ r f ,

0≤x≤L

DΔCH , t

(

(12)

Several approaches have been undertaken to express the growth rate, G, as a function of the supersaturation.25−28 Here, the first derivative of the desupersaturation curve is used for the calculation of ki. Hence, the deposition rates (mass transfer due to the crystallization) vary for different cooling temperatures and flow rates, and the desupersaturation profiles are specific for a particular condition.29 Therefore, G, ki, and the crystal growth parameter should be calculated individually for each combination of temperature and flow rate. 2.4. Heat Transfer. Similar to the mass transfer equation, the heat transfer governing equation with the initial and boundary conditions is defined by eq 13.

Keq is the corresponding equilibrium distribution coefficient of the local concentration of impurities at a given time (relative concentration of impurities, RelCimpt), msolinc,t is the mass of the included solution at a given time, ms,t is the mass of deposited solids in the crystal layer at a given time, Cimp,t is the local concentration of impurities at a given time, and KG represents the entrapment of impurities due to the crystal growth (eq 9). Jiang et al.3,8 proposed an effective distribution coefficient for entrapment of impurities due to crystal growth in melt crystallization. The same concept can be employed for falling film solution crystallization, in which the rapid crystal growth rate causes the formation of an imperfect crystal layer (porous). KG = 1 −

drs k (C − Csat )n 1 dm = = i r dt ρs A dt ρs

(11) 5022


Article

Industrial & Engineering Chemistry Research 6. Using the second and third order Runge Kutta scheme to evaluate the expected values of the temperature and concentration. 7. Updating the feed tank concentration for recirculation. 8. Calculating the new solid thickness over the length of the column for heat transfer. 9. Returning to step 3 for the next time step until the total simulation time was reached. 10. Calculating yield and purity over time. Figure 3 shows the numerical algorithm used to solve the model.

Figure 4. Change in average purity for system I at different flow rates and cooling temperatures.

Figure 3. Numerical algorithm for solving the coupled partial differential equations for the falling film crystallizer with a number of columns in series.

Three different solution systems, which were prepared for the experiments,2 were used for this modeling study, which are listed in Table 1. The physical properties, initial, and boundary conditions are listed in the Supporting Information.

Figure 5. Change in average purity for system II at different flow rates and cooling temperatures.

Table 1. Crystallization Solution Systems for the Falling Film Experiments systems system I

main compound (initial concentration, w %)

impurity (initial concentration, w %)

system II

acetaminophen (95%) fenofibrate (98%)

metacetamol (5%) fenofibric acid (2%)

system III

fenofibrate (98%)

fenofibric acid (2%)

solvent ethanol and water (50:50 v %) ethyl acetate and ethanol (70:30 v %) ethanol (100%)

4. MODELING RESULTS AND MODEL VALIDATION BY THE EXPERIMENTAL RESULTS 4.1. Average Purity of the Deposited Crystal Layer. The two key performance indicators of the process are final purity of the crystal layer and yield. Figures 4, 5, and 6 show the modeling results for average purity of the deposited layer for the systems at different flow rates and cooling temperatures. It is obvious for all three systems that the highest flow rate improved the purity. This improvement in purity is the result of

Figure 6. Change in average purity for system III at different flow rates and cooling temperatures.

5023


Article

Industrial & Engineering Chemistry Research a thinner boundary layer at the higher flow rate because the concentration of the impurities (the concentration profile of the impurities in the radial direction) is lower in the thinner boundary layer on the growing face of the crystal. Therefore, the rejection rate is higher, and fewer impurities are entrapped in the forming crystal lattice. Conversely, the lowest level of purity corresponds to the lowest flow rate of 5 mL·min−1, when the rejection rate of impurities on the boundary layer is lower. For system I (Figure 4) with a starting amount of 5% impurity content at 0 °C, flow rates of 5, 20, and 30 mL·min−1 yielded average final purities of 96.6%, 96.8%, and 97.1%, respectively. The initial nucleation phase periods had no meaningful level of purity for a yet (not) deposited layer. A comparison between systems II and III shows the effect of solvent selection, in which the growth rate (and effective distribution coefficient) can be manipulated using a variety of solvents and mixtures. The results for the different temperatures and flow rates show no major effect on the final purity, although the final purity is significantly improved compared to the initial amount of impurities in the solution. The solvent selection and crystallization solution system design have a great effect on the final purity, due to alterations of the crystal growth rate, in which a higher growth rate will produce a higher impurity inclusion rate (i.e., comparing systems II and III). 4.2. Relative Concentration of Impurity and Effective Distribution Coefficient. The falling film crystallization experiments and this modeling study were performed in the recirculation mode. The concentration of the feed tank is changed due to the crystallization rate and the flow rate of the feed. The crystal growth rates and mass transfer coefficients were calculated on the basis of this desupersaturation profile. Because the concentration of solution that enters the crystallizer column changes over time, the crystal growth parameters were not constant and were dynamically estimated for each time scale. Because this cooling crystallization does not have a 100% purification yield, the concentration of impurity is changing over process time, with impurity inclusion in the crystal layer. However, the proportional rate of change of concentration of API to impurity is not constant. The modeling results for the change of concentration of APIs and impurities in the feed tank for all systems can be found in the Supporting Information. The relative concentration of impurities to the concentration of the API had a significant effect on the distribution coefficient and overall inclusion of impurities in the crystals. Figures 7, 8, and 9 present the profiles of the relative concentration of impurities for the systems with different flow rates and cooling temperatures. The initial portion of the profiles shows a large slope, where the deposition rate of the main API was high due to the higher supersaturation degree of the API at the beginning and the rapid deposition of crystals, which resulted in a rapid decrease in the concentration of the API. However, the concentration of impurities did not change at the same rate. A turning point occurs at the later stage when the crystal growth rate and deposition rate decrease and a high concentration of impurities is included. For instance, for system I, at a flow rate of 20 mL·min−1 and a cooling temperature of 0 °C, the initial relative concentration of impurities is 5% wimp·wAPI−1. The relative concentration increases to 8.3% wimp·wAPI−1 and then decreases to a steady level. At a high relative concentration of impurities, the amount of impurity inclusion changes due to an increase in the equilibrium distribution coefficient and the higher concen-

Figure 7. Relative concentration of impurities for system I at different flow rates and cooling temperatures.

Figure 8. Relative concentration of impurities for system II at different flow rates and cooling temperatures.

Figure 9. Relative concentration of impurities for system III at different flow rates and cooling temperatures.

5024


Article

Industrial & Engineering Chemistry Research tration of solution entrapped in crystal cavities and between crystal cleavages (solvent inclusion). The results demonstrate the need to remove excess amounts of impurities during the crystallization process to restrict the maximum relative concentration of impurities and achieve better purification. Distribution coefficients are the key factors for purification. The effective distribution coefficient is a kinetic parameter, and it depends on many processing parameters as well as the equipment size, design, and geometry. Figures 10, 11, and 12

Figure 12. Distribution coefficients for system III at different flow rates and cooling temperatures.

transfer and growth parameters) on the effective distribution coefficient are significant. For instance, Figure 10 shows the distribution coefficients for system I, in which the bottom curve represents the equilibrium condition. The value for the equilibrium distribution coefficient for 5% wimp·wAPI−1 is 0.14, and for the lowest effective distribution coefficient (20 mL· min−1, 0 °C), the value is 0.52, which is approximately 3.7 times greater; the ratio is 4.9 for the top curve (5 mL·min−1, 0 °C). Better purification performance is expected for the system and process conditions with lower effective distribution coefficients. In the case of system II (Figure 11), the effective distribution coefficient for the experiment conducted at a flow rate of 10 mL·min−1 and a cooling temperature of 10 °C was 0.85, which is 31.5 times higher than the equilibrium distribution coefficient; therefore, the purification was not satisfactory in this case. This high effective distribution coefficient is due to the very high growth rate of the API crystal under those system and process conditions, which causes the entrapment of many impurities in the crystal matrix and deposited layer. 4.3. Yield of the Process. Yield is the other main performance indicator of the process. The modeling results for yield of the falling film crystallization process for the systems at different flow rates and cooling temperatures is presented in Figures 13, 14, and 15. The overall yield for all of the systems is approximately 70%. This relatively low yield is due to the recirculation of the feed and maintaining the feed temperature at 60 °C. The deposited crystal layer causes a resistance to heat transfer, which results in a higher liquid film temperature toward the end of the experiments. The modeling results for the layer thickness can be found in the Supporting Information. The heat transfer model results in the following section demonstrate this effect. Therefore, toward the end of experiments, the concentration is much lower than the starting concentration, and the equilibrium temperature of the falling liquid film is also higher. These effects lead to a significantly lower supersaturation ratio, which is the driving force for the cooling crystallization. Therefore, the remaining dissolved APIs can not be deposited in the crystal layer, exit the crystallizer in a cool and low supersaturation state, enter the feed tank, and are heated up and recirculated in the crystallizer again. An effective and practical solution to enhance the yield and continue harvesting crystals would be to lower the feed tank temperature during the middle of the experiments.

Figure 10. Distribution coefficients for system I at different flow rates and cooling temperatures.

Figure 11. Distribution coefficients for system II at different flow rates and cooling temperatures.

show the effective and equilibrium distribution coefficients for the systems at different flow rates and cooling temperatures. The curves are calculated on the basis of actual experimental results, which are included in the graphs, in which the final purity of the harvested crystals were measured by the HPLC technique, and the graphs were generated on the basis of the initial concentration of impurities. The experimental values are average of three repeats with significant error of less than 5%(error bars are not shown here). The trend of distribution coefficients from experimental values was used to fit the curves for the conditions where experimental values for intermediate impurity concentrations were not available. For instance, for Figure 10, the detail results for flow rate of 5 mL/min and 10 °C was used to fit a curve based on the proportional difference of the 5% impurity concentration for flow rate of 20 mL/min and 10 °C. The effects of temperature and flow rate (mass 5025


Article


Figure 13. Overall yield profile for system I at different flow rates and cooling temperatures.

Figure 15. Overall yield profile for system III at different flow rates and cooling temperatures.

Figure 14. Overall yield profile for system II at different flow rates and cooling temperatures.

Figure 16. Axial temperature profile of system I at a flow rate of 20 mL·min−1 and a cooling temperature of 0 °C at different process times.

4.4. Axial Temperature Profile. The axial temperature profiles for systems I and II at a flow rate of 20 mL·min−1 and a cooling temperature of 0 °C at different processing times are shown in Figures 16 and 17. At the start of the process, the feed temperature is 60 °C before introducing the feed to the column. The top line in the graphs is a hypothetical temperature profile for time 0 to show the starting point of the temperature. By the start of the experiment, when the hot solution is uniformly distributed over the column and begins to descend, heat transfer occurs, and the liquid film cools to an equilibrium temperature depending on the cooling liquid inside the pipe. The body of pipe is metal and has high conductivity and low heat transfer resistance; therefore, the temperature of the liquid film approaches 0 °C (in this case) at the top of the column. By developing a less conductive crystal layer, the heat transfer rate decreases, and the equilibrium temperature shifts to the lower axial position but eventually reaches equilibrium and a steady state condition. The final temperature depends on the thickness of the deposited layer. For instance, for system II, in which the final average thickness of the layer was almost

double that of the final average thickness of system I (5.5 vs 2.9 mm), the final equilibrium temperature was higher (5.14 vs 3.2 °C), although the ratio was not linear. 4.5. Performance of a Multicolumn System. A simple approach to improving the purity of the product would be to employ a multicolumn setup in series. In this configuration, a solid crystal layer deposited on the column could perform secondary purification operations. Clean solvent could be introduced to dissolve or partially dissolve the layer to deliver purified product for the next operation, producing another falling film stage for further purification. The experimental and simulation results showed that single stage purification significantly decreased the level of impurities, for instance in the case of system I at a flow rate of 20 mL·min−1 and a cooling temperature of 0 °C, the initial concentration of impurities decreased from 5% wimp·wAPI−1 to 3.2% wimp·wAPI−1. If this layer were dissolved in a fresh solvent, a solution with a 3.2% wimp· wAPI−1 impurity concentration could be further purified in a consecutive stage(s). An experiment and simulation study was 5026


Article


deposited layer from the actual experiment at the end of the crystallization time (12 600 s) was the same as that obtained in the simulation. The average purities of the deposited layers at the different time intervals were lower than those obtained in the simulation results. This variation could be due to the effective distribution coefficient calculation, which was based on the final purity and final amount of impurity inclusion, although the impurity inclusion level and mechanism during layer deposition over time are not fully understood. However, the overall simulation results are consistent and reliably confirm the experimental values. 4.6. Variable Temperature Control for Improving Yield and Final Purity. As mentioned in Section 4.3, the average process yield is around 70%. Section 4.4 hypothesized the effect of liquid film temperature at a later stage of experiments on this relatively low yield. The final purity results (Section 4.2) show a considerable difference between achievable purity from the process compared with the theoretical purity, based on equilibrium distribution coefficients. The effect of variable feed and column temperatures on the yields and purities were studies for system I here (acetaminophen (95%) and metacetamol (5%) in an ethanol/ water mixture, flow rate of 20 mL·min−1). Two scenarios were tested: (1) constant column temperature (0 °C) and variable feed temperature from 60 to 20 °C (ramp of −0.5 °C/min); (2) variable column temperature from 30 to 0 °C (ramp of −0.5 °C/min) and variable feed temperature from 60 to 30 °C (ramp of −0.5 °C/min). For scenario 1 (variable feed temperature), the liquid falling film enters the column at low temperature toward the end of the process, when a less conductive crystal layer on the column lowers the heat transfer rate. Therefore, the lower heat transfer rate is sufficient to effectively cool down the liquid film from the top of the column and provides higher supersaturation and consecutively contributes in further crystallization. The yield improved from 69% to 79%, but the final purity did not change significantly. As explained in the previous sections, the final purity is more dependent on the growth rate, and the high growth rate due to the large temperatures difference for most of the process time causes most of the impurity inclusion in the crystal layer. For scenario 2 (variable feed and column temperatures), the degree of supersaturation was almost constantly maintained for most of the processing time by keeping a temperature difference of 30 °C (feed and column temperatures), although the solubility for the range of temperatures is not linear and so the degree of supersaturation. The yield improved from 69% to 78% because of the aforementioned temperature effect. The controlled growth rate, which was manipulated by temperature, had significant effect on improving the final purity from 96.8% to 97.9%. Table 2 shows the comparison of final purities and yields for fixed and variable temperatures of feed and column. Combinations of variable cooling temperature and processing time (ramp of the change of temperatures) can be tried in

Figure 17. Axial temperature profile of system II at a flow rate of 20 mL·min−1 and a cooling temperature of 0 °C at different process times.

performed for a solution with the same API concentration (acetaminophen in this case) and a 3% wimp·wAPI−1 impurity concentration at the second stage of falling film crystallization. The simulation and experimental results are shown in Figure 18. The fitted effective distribution coefficient for the system

Figure 18. Average purity of system I at a flow rate of 20 mL·min−1 and a cooling temperature of 0 °C for a multistage process (stage 1 with 5% initial impurity, stage 2 with 3% initial impurity).

with a 3% impurity concentration (Figure 10) was used for the simulation. The crystallization experiment was repeated with different processing time intervals, and samples from the deposited layer were tested by HPLC. The final purity of the

Table 2. Comparison of Final Purities and Yields for Fixed and Variable Temperatures of Feed and Columna

a

scenario

feed temperature

column temperature

yield (%)

final purity (w %)

fixed temperature gradual feed cooling (column fix) gradual feed and column cooling

60 °C 60−20 °C (−0.5 °C/min) 60−30 °C (−0.5 °C/min)

0 °C 0 °C 30−0 °C (−0.5 °C/min)

69 79 78

96.8 97.0 97.9

System I, acetaminophen (95%) and metacetamol (5%) in an ethanol/water mixture, flow rate of 20 mL·min−1. 5027



■

order to improve yield and final purity. This strategy can be extremely helpful for fast growing systems, such as system II, where manipulation of the high growth rate can reduce dendritic crystal formation and lower the amount of impurity inclusion. For experimental works, the high initial column temperature can be challenging for the nucleation phase. In this case, after the first stage of nucleation phase, the temperatures can be tuned to control the crystallization rate.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We would like to acknowledge Novartis International AG for its generous funding of this research.

5. CONCLUSIONS

■

A numerical study of the falling film crystallization process over a cold column in recirculation mode was conducted to predict and model the solute concentration, temperature profile, crystal layer thickness, yield of the process, and final purity of the product as a function of time. The mathematical model included the heat and mass transfer, crystal growth rate, and impurity inclusion. The simulation results were compared with experimental values to evaluate the predictive capability of the mathematical model. A range of flow rates and cooling temperatures for three different solution systems was used for the model input. The mathematical model accurately predicted the system performance. In addition, the model was used to predict a multistage process for further purification of an API. The mathematical model can be used for further optimization of the process and to determine the results for a combination of process parameters. Although the physical limitations of the crystallization process, such as the nucleation phase, restrict the extrapolation of results from the model, combinations of the processing parameters within the design space may be reliable. The results are also based on the proposed crystallizer design, size, and geometry. In the case of scaling up (or down) or a change in geometry, new experiments should be conducted to regress the model inputs and validate the model results. A comparison between systems II and III showed the effect of the growth rate (and distribution coefficient) on the overall performance and, in particular, the purification of the system, which can be manipulated by varying the solvents and mixtures. In general, the results for the different temperatures and flow rates showed a major effect on the final purity due to the mass transfer and growth rate terms. The solvent selection and crystallization solution system design had a greater effect on the final purity, due to the alteration of the crystal growth rate, in which a higher growth rate produced a higher impurity inclusion rate (i.e., comparing systems II and III). On the basis of the simulation results, optimization of the system for parameters such as the overall yield can also be achieved by dynamically changing the input parameters, such as the temperature of the feed tank during an experiment.

■

Article

NOMENCLATURE A Area (m2) C Concentration (kg·kg−1) Cp Heat capacity (j·kg−1·K−1) D Diffusivity (m2·s−1) F Flow rate (m3·s−1) G Crystal growth rate (m·s−1) H Heat transfer coefficient (W·m−2·K−1) J Mass flux (kg·m−2·s−1) K Distribution coefficient (−) k Thermal conductivity (W·m−1·K−1) kd Mass transfer coefficient (m·s−1) ki Integration growth rate (m·s−1) L Length of column (m) M Mass (kg) n Growth rate parameter (kinetic exponent) (−) RelCimp Relative concentration of impurities (−) r Radius (m) sh Sherwood number from Ranz-Marshall equation (−) Sj Number of stages (column) (−) T Temperature (K) t Time (s) u Velocity (m·s−1) U Overall heat transfer coefficient (W·m−2·K−1) V Volume (m3) x Axial position (m) Greek Letters

α μ ρ ϕ

Thermal diffusivity (m2·s−1) Dynamic viscosity (Pa·s) Density (kg·m−3) Column diameter (m)

Subscripts

0 col eq f G H imp m s sat sol solinc solv tank t w

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b00057. Initial values and boundary conditions; experimental results; solubility curves; modeling results for desupersaturation profile of APIs in the feed tank; modeling results for changes in the concentration of impurities in the feed tank; modeling results for change in deposited layer thickness over time. (PDF) 5028

Initial state Column Equilibrium Fluid Growth Host compound Impurity Metal Solid, Solute Saturation Solution Solvent inclusion Solvent Feed tank Time Cooling water DOI: 10.1021/acs.iecr.6b00057 Ind. Eng. Chem. Res. 2016, 55, 5019−5029

Article


■

Mechanisms and Kinetics. Theor. Found. Chem. Eng. 2003, 37, 137− 143. (22) Zhou, L.; Su, M.; Benyahia, B.; Singh, A.; Barton, P. I.; Trout, B. L.; Myerson, A. S.; Braatz, R. D. Mathematical modeling and design of layer crystallization in a concentric annulus with and without recirculation. AIChE J. 2013, 59, 1308−1321. (23) Zhang, H.; Quon, J.; Alvarez, A. J.; Evans, J.; Myerson, A. S.; Trout, B. Development of Continuous Anti-Solvent/Cooling Crystallization Process using Cascaded Mixed Suspension, Mixed Product Removal Crystallizers. Org. Process Res. Dev. 2012, 16, 915−924. (24) Ulrich, J.; Neumann, M. Purification by solid layer melt crystallization. J. Therm. Anal. 1997, 48, 527−533. (25) Cornel, J.; Mazzotti, M. Estimating Crystal Growth Rates Using in situ ATR-FTIR and Raman Spectroscopy in a Calibration-Free Manner. Ind. Eng. Chem. Res. 2009, 48, 10740−10745. (26) Mitchell, N. A.; Ó ’Ciardhá, C. T.; Frawley, P. J. Estimation of the growth kinetics for the cooling crystallisation of paracetamol and ethanol solutions. J. Cryst. Growth 2011, 328, 39−49. (27) Qiu, Y.; Rasmuson, Å. C. Crystal growth rate parameters from isothermal desupersaturation experiments. Chem. Eng. Sci. 1991, 46, 1659−1667. (28) Zipp, G. L.; Rodriguez-Hornedo, N. Determination of crystal growth kinetics from desupersaturation measurements. Int. J. Pharm. 1989, 51, 147−156. (29) 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, 1625−1628.

REFERENCES

(1) Wong, S. Y.; Chen, J.; Forte, L. E.; Myerson, A. S. Compact Crystallization, Filtration, and Drying for the Production of Active Pharmaceutical Ingredients. Org. Process Res. Dev. 2013, 17, 684−692. (2) Yazdanpanah, N.; Ferguson, S. T.; Myerson, A. S.; Trout, B. L. Novel Technique for Filtration Avoidance in Continuous Crystallization. Cryst. Growth Des. 2016, 16, 285−296. (3) Jiang, X.; Hou, B.; He, G.; Wang, J. Falling film melt crystallization (I): Model development, experimental validation of crystal layer growth and impurity distribution process. Chem. Eng. Sci. 2012, 84, 120−133. (4) Jiang, X.; Li, M.; He, G.; Wang, J. Research Progress and Model Development of Crystal Layer Growth and Impurity Distribution in Layer Melt Crystallization: A Review. Ind. Eng. Chem. Res. 2014, 53, 13211−13227. (5) Kim, K.-J.; Ulrich, J. A study on effective thermal conductivity of crystalline layers in layer melt crystallization. J. Phys. D: Appl. Phys. 2002, 35, 1080. (6) Beierling, T.; Gorny, R.; Sadowski, G. Modeling Growth Rates in Static Layer Melt Crystallization. Cryst. Growth Des. 2013, 13, 5229− 5240. (7) Beierling, T.; Osiander, J.; Sadowski, G. Melt crystallization of isomeric long-chain aldehydes from hydroformylation. Sep. Purif. Technol. 2013, 118, 13−24. (8) Jiang, X.; Xiao, W.; He, G. Falling film melt crystallization (III): Model development, separation effect compared to static melt crystallization and process optimization. Chem. Eng. Sci. 2014, 117, 198−209. (9) LU, D.-P.; LIU, Y.; HE, G.; JIANG, X.-B. Fractal Analysis for Crystal Layer Structure and Impurity Distribution Research of Melt Crystallization. Fractals 2015, 23, 1540002. (10) Beierling, T.; Micovic, J.; Lutze, P.; Sadowski, G. Using complex layer melt crystallization models for the optimization of hybrid distillation/melt crystallization processes. Chem. Eng. Process. 2014, 85, 10−23. (11) Chianese, A.; Santilli, N. Modelling of the solid layer growth from melt crystallizationthe integral formulation approach. Chem. Eng. Sci. 1998, 53, 107−111. (12) Guardani, R.; Belline, A. Application of mathematical modelling in a two-component layer crystallisation process. Chem. Eng. Technol. 1997, 20, 495−501. (13) Guardani, R.; Neiro, S. M. S.; Bülau, H.; Ulrich, J. Experimental comparison and simulation of static and dynamic solid layer melt crystallization. Chem. Eng. Sci. 2001, 56, 2371−2379. (14) Chen, X. D.; Wu, W. D.; Chen, P. An analytical relationship of concentration-dependent interfacial solute distribution coefficient for aqueous layer freeze concentration. AIChE J. 2015, 61, 1334−1344. (15) Kim, S.; Seo, M.; Tak, M.; Kim, W.; Yang, D.; Kang, J. Effects of sweating time and cooling strategy on purification of N-vinyl-2pyrrolidinone using a melt crystallizer. Korean J. Chem. Eng. 2013, 30, 1997−2000. (16) Kim, K. J. Purification of Phosphoric Acid from Waste Acid Etchant using Layer Melt Crystallization. Chem. Eng. Technol. 2006, 29, 271−276. (17) Kim, K.-J.; Ulrich, J. Impurity Distribution in a Solid−Liquid Interface during Static Layer Crystallization. J. Colloid Interface Sci. 2002, 252, 161−168. (18) Taran, A. L.; Nosov, G. A.; Myasnikov, S. K.; Kholin, A. Y. Kinetics of Crystallization of Binary Melts of Eutectic-Forming Substances. Theor. Found. Chem. Eng. 2004, 38, 164−168. (19) Myasnikov, S. K.; Uteshinsky, A. D. Melt entrapment during growth of crystal layer with nonplanar interface and diffusive transport of captured impurities out of layer pores. J. Cryst. Growth 2005, 275, e39−e45. (20) Parisi, M.; Chianese, A. The crystal layer growth from a wellmixed melt. Chem. Eng. Sci. 2001, 56, 4245−4256. (21) Myasnikov, S. K. Transport of Impurities out of a Two-Phase Crystal Layer into a Melt under the Effect of a Temperature Gradient: 5029


Mathematical Modeling of Layer Crystallization on a Cold Column

Recommend Documents