Crystallization of Benzoic Acid in a Mixed Suspension, Mixed Product Removal Crystallizer

Crystallization is a key separation process in chemical engineering. Hundreds of millions of tons of materials are produced by crystallization yearly. Some examples of materials produced by crystallization include sodium chloride, sodium sulphates and sucrose. A main motivation for using crystallization is that separating a product from a solution by crystallizing it out of solution usually requires less energy than evaporating the solvent.

Crystallization takes place in two steps. First, a crystal nuclei must form (nucleation). Second, the nuclei grows by deposition of more crystal material (growth). The driving force for crystallization is the degree of supersaturation. Supersaturation can be achieved for example by cooling the solution. Upon cooling, the solution becomes supersaturated, and crystallization can begin. The key objective in designing and operating a crystallizer is to control the crystal size distribution. For example, in the crystallization of a pharmaceutical product, the size of the crystals formed will determine the time the product will take to dissolve in the stomach of a patient.

In this example, the crystallization of benzoic acid is modeled in an idealized continuous crystallizer known as a “mixed suspension, mixed product removal reactor”. A 0D model is defined, and a population balance equation is solved to calculate the crystal size distribution and the mean particle size under different conditions. The model is based on a paper by Morris and others (Ref. 1). Figure 1 shows a sketch of the crystallizer setup.

Figure 1: Sketch of the MSMPR setup. The crystallizer is kept at a temperature below that of the feed to create supersaturation.

Model Definition

The crystallization of benzoic acid is modeled in 0D. The crystallizer and the feed vessel are both assumed to be perfectly mixed and operated in steady state. Kinetic parameters fitted to experimental data in Ref. 1 are used.

For an ideal MSMPR crystallizer, the analytical solution to the population balance equation is

(1)

where n is the population density, n0 is the nuclei population density, L is the crystal size, G is the crystal growth rate, and τ is the mean residence time. The mass-based mean crystal size distribution is found by taking the ratio of the fourth and third moments of n as defined by Equation 1. The nuclei population density is obtained from

(2)

where B0 is the nucleation rate. The nucleation rate is given by

(3)

where kb is a nucleation rate coefficient, MT is the MSMPR suspension crystal weight fraction, ΔC is the driving force for crystallization, and j and b are dimensionless parameters. The driving force for crystallization is given by

(4)

where C is the benzoic acid mass fraction (g benzoic acid/kg solvent) in the MSMPR crystallizer, and C*(T) is the saturation mass fraction.

The saturation mass fraction is obtained from a polynomial expression determined experimentally for a mixture of benzoic acid in 1.5:1 (w/w) water/ethanol:

(5)

A graphical representation of C* is found in Figure 2.

Figure 2: The saturation mass fraction of benzoic acid as a function of temperature.

The growth rate is obtained from

(6)

where kg is the crystal growth rate coefficient and g is a fitted kinetic parameter. In turn, the crystal growth rate coefficient is obtained from

(7)

where kg0 is the Arrhenius expression pre-exponential factor, and Eg is the activation energy for crystal growth.

Solving for the crystal concentration allows the population density to be determined. A mass balance on the crystallizer gives the equation

(8)

Here, C0 is the concentration of benzoic acid in the feed. Furthermore, the suspension crystal weight fraction is related to the nuclei population density and the growth rate through

(9)

where kv and ρc are the volumetric shape factor and density of the crystals, respectively.

The crystal concentration is obtained by solving for the difference of these two definitions of MT, that is

(10)

Equation 10 is solved for in a node.

The steady-state mass-weighted mean crystal size, L4,3, is calculated from the ratio of the fourth and third moments of the crystal size distribution, using the following expression:

(11)

Results and Discussion

Crystal size distributions are determined for the experimental setup giving the largest and smallest crystal sizes, and for sweeps from low to high reactor temperature, and low to high residence time. Finally, the effect on the mean particle size of both reactor temperature, and of residence time, is examined.

Crystal size distributions

Figure 3: Crystal size distributions for the experimental conditions where the smallest and largest particles are produced.

Figure 3 shows the crystal size distributions for the combinations of temperature and residence time that results in production of the smallest and the largest particles. These results agree well with those from Ref. 1.

Figure 4: Crystal size distributions for increasing crystallizer temperatures.

Figure 2 shows that the saturation mass fraction of benzoic acid increases with temperature. This reduces the driving force for crystallization, ΔC. At the same time, the growth rate of the crystals increases with increasing temperature, according to Equation 6. The effect on the crystal size distribution is that the crystals that do form, grow to larger average sizes. Increasing the residence time also has a similar effect, as seen in Figure 5.

Figure 5: The effect of residence time on the crystal size distribution.

Mean particle size

Figure 6 and Figure 7 show the average crystal size for different temperatures and residence times, respectively.

Figure 6: The effect of temperature on the average particle size. The residence time is 45 min.

Figure 7: The effect of residence time on the average particle size. The crystallizer temperature is 30°C.

Notes About the COMSOL Implementation

The original paper (Ref. 1) defined crystal concentrations in both g/mL solvent and g/g solvent, but for simplicity, this model defines concentrations only in g/kg solvent. This does not affect the resulting crystal size distributions.

References

1. G. Morris, G. Power, S. Ferguson, M. Barrett, G.Hou, and B. Glennon, “Estimation of Nucleation and Growth Kinetics of Benzoic Acid by Population Balance Modeling of a Continuous Cooling Mixed Suspension, Mixed Product Removal Crystallizer,” Org. Process Res. Dev., vol. 19, pp. 1891–1902, 2015.

2. J.F. Richardson, J.H. Harker, and J.R. Backhurst, Coulson & Richardson’s Chemical Engineering, vol. 2, 5th edition, Elsevier, 2002.

Chemical_Reaction_Engineering_Module/Mixing_and_Separation/benzoic_acid_crystallization

Modeling Instructions

A 0D model using a Global Equations interface will be used to describe the crystallization process.

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Import the parameters for the model.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file benzoic_acid_crystallization_parameters.txt.

In the toolbar, click

Define two parameter cases to simulate the conditions that give rise to the smallest and the largest crystals.

In the window for , type Experiment Giving Largest Crystals in the text field.

In the toolbar, click

In the window for , type Experiment Giving Smallest Crystals in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Description
T	0[degC]	MSMPR operating temperature
tau	20 [min]	Residence time

Define a Function representing the saturation mass fraction.

Saturation Mass Fraction

In the toolbar, click

and choose .

In the window for , type Saturation Mass Fraction in the text field.

In the text field, type C_star.

Locate the section. In the text field, type 2.03e-5*T^4 + 2.97e-4*T^3 + 4.70e-2*T^2 + 1.43*T + 24.71.

In the text field, type T.

Locate the section. In the text field, type g/kg.

In the table, enter the following settings:


Argument	Unit
T	degC

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
T	-15[degC]	40[degC]	K

Click

The saturation mass fraction increases with temperature. Decreasing the temperature is a key driving force for crystallization.

To simplify the model setup, import the required variable expressions from a separate file.

Definitions

Variables 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file benzoic_acid_crystallization_variables.txt.

Define the global equation that will be solved to obtain the benzoic acid mass fraction.

Global ODEs and DAEs (ge)

Global Equations 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	f(u,ut,utt,t) (1)	Initial value (u_0) (1)	Initial value (u_t0) (1/s)	Description
C	M_T_MB-M_T_MR	C_star(T)+1[g/kg]	0	Benzoic acid mass fraction

Locate the section. Click

In the table, enter the following settings:


Dependent variable quantity	Unit
Custom unit	g/kg

Click

In the table, enter the following settings:


Source term quantity	Unit
Custom unit	g/kg

Study 1: Extreme Cases

In the window, click .

In the window for , type Study 1: Extreme Cases in the text field.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Click

In the toolbar, click

Plot the crystal size distributions for the experiments giving the smallest and the largest crystals .

Results

Grid 1D: Extreme Cases

In the toolbar, click

and choose .

In the window for , type Grid 1D: Extreme Cases in the text field.

Locate the section. In the text field, type L.

Locate the section. From the list, choose .

Locate the section. In the text field, type L_max.

CSD for Extreme Cases

In the toolbar, click

In the window for , type CSD for Extreme Cases in the text field.

Locate the section. From the list, choose .

Click to expand the section. From the list, choose .

In the text area, type Volume-based Mean Crystal Size Distribution.

Locate the section.

Select the check box. In the associated text field, type Crystal Size Distribution (%).

Locate the section. Select the check box.

Locate the section. From the list, choose .

Line Graph 1

In the toolbar, click

In the window for , locate the section.

In the text field, type (n0 * exp(- L / (G*tau))*L^4)/ integrate((n0 * exp(- L / (G*tau)))*L^3, L, 0, L_max).

From the list, choose .

Locate the section. From the list, choose .

In the text field, type L.

From the list, choose .

Select the check box. In the associated text field, type Crystal size.

Click to expand the section. Select the check box.

From the list, choose .

In the table, enter the following settings:


Legends
T =20\deg C, \tau = 45 min
T = 0\deg C, \tau = 15 min

In the toolbar, click

CSD for Extreme Cases

In the window, click .

In the window for , locate the section.

Select the check box.

In the text field, type 10.

In the toolbar, click

Next, we’ll set up a study to calculate the effect of temperature on the crystal size. For this study, the longest residence time will be used.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 2: Temperature Sweep

In the window, click .

In the window for , type Study 2: Temperature Sweep in the text field.

Step 1: Stationary

In the window, under click .

In the window for , click to expand the section.

Select the check box.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
T (MSMPR operating temperature)	range(0,1,30)	degC

Solution 5 (sol5)

In the toolbar, click

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type 0.0001.

Reduce the solver tolerance to avoid numerical noise in the temperature sweep results.

In the toolbar, click

Results

Grid 1D: Temperature Sweep

In the window, right-click and choose .

In the window for , type Grid 1D: Temperature Sweep in the text field.

Locate the section. From the list, choose .

CSD for Temperature Sweep

In the window, right-click and choose .

In the window for , type CSD for Temperature Sweep in the text field.

Locate the section. From the list, choose .

From the list, choose .

In the list, choose , , and .

Line Graph 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

CSD for Temperature Sweep

In the window, click .

In the toolbar, click

Next, we’ll set up a study to calculate the effect of residence time on the crystal size. For this study, the highest temperature will be used.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 3: Residence Time Sweep

In the window, click .

In the window for , type Study 3: Residence Time Sweep in the text field.

Step 1: Stationary

In the window, under click .

In the window for , locate the section.

Select the check box.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
tau (Residence time)	range(15,1,45)	min

In the toolbar, click

Results

Grid 1D: Residence Time Sweep

In the window, right-click and choose .

In the window for , type Grid 1D: Residence Time Sweep in the text field.

Locate the section. From the list, choose .

CSD for Residence Time Sweep

In the window, right-click and choose .

In the window for , type CSD for Residence Time Sweep in the text field.

Locate the section. From the list, choose .

In the toolbar, click

Mean Particle Size for Different Operating Temperatures

In the toolbar, click

and choose .

Finally, create two plots showing the effect on average particle size of temperature, and residence time, respectively.

In the window for , type Mean Particle Size for Different Operating Temperatures in the text field.

Locate the section. From the list, choose .

Global 1

In the toolbar, click

The mean particle size plots do not use a Grid dataset. For this reason, L is not defined. We will instead use the formal integration constant parameter L_int.

In the window for , locate the section.

In the table, enter the following settings:


Expression	Unit	Description
integrate(n0 * exp(- L_int / (Gtau))L_int^4, L_int, 0, L_max) / integrate((n0 * exp(- L_int / (Gtau)))L_int^3, L_int, 0, L_max)	um	Mean Particle Size

Locate the section. From the list, choose .

In the text field, type T.

From the list, choose .

In the toolbar, click

Mean Particle Size for Different Residence Times

In the window, right-click and choose .

In the window for , type Mean Particle Size for Different Residence Times in the text field.

Locate the section. From the list, choose .

Global 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type tau.

From the list, choose .

In the toolbar, click