ABSTRACT

Mefenamic acid (MA) is a high-dose, anti-inflammatory, analgesic agent that is widely prescribed for pain related to menstrual disorders. It has some negative properties, such as a high hydrophobicity with a propensity to stick to surfaces, and possess great problems during granulation and tableting. Crystallization kinetics was investigated for mefenamic acid. Availability of data on the kinetics of crystal growth is very important for the development and operation of industrial crystallisation processes. The experiments for the measurement of crystal growth kinetics were carried out using the desupersaturation curve technique based on the measurement of the solution concentration versus time in a seeded isothermal batch experiment. To predict the optimum parameters (b, k_b, g, k_g) for the nucleation and growth kinetics from the desupersaturation curve obtained, the Population Balance Modelling was used and solved by the method of moments. The initial values for the optimisation problem were estimated by using the approach developed by Garside et al. (1982) [1].

1. Introduction

Mefenamic acid [N-(2,3-xyly) anthranilic acid] is a potent prostaglandin synthetase inhibitor that is widely used as a non-steroidal anti-inflammatory and analgesic-antipyretic drug. Mefenamic acid (MA) is used for the relief of postoperative traumatic inflammation and swelling, antiphlogistic analgesic treatment of rheumatoid arthritis, for antipyresis in acute respiratory tract infection and arthritis. Mefenamic acid is available in the form of tablets, capsules, and a pediatric suspension. Mefenamic acid is a problematic drug in granulation, tableting, and dissolution due to its poor solubility, hydrophobicity, and tendency to stick to surfaces.

It has a prevalent usage in drug formulation. It is a poorly soluble drug in an aqueous media, sticking to any type of surface, and is not easily handled in the granulation and tableting processes. As is the case for many poorly soluble drugs, dissolution of mefenamic acid is a problem. To facilitate tableting, enhance dissolution rates, and develop a stable and reproducible dosage form, investigation of the physicochemical properties of mefenamic acid is necessary.

Pharmaceutical drugs are commonly crystalline materials and are therefore subject to polymorphism. As a result of different arrangement and or conformation, the polymorphic solids have different unit cells and hence exhibit different physical properties, including intermolecular interaction, crystalline disorder, a change in free energy, thermodynamic activity, solubility, dissolution rate, stability, and bioavailability.

Many of these properties, in particular stability and solubility, have an impact on the development potential of a drug molecule and it is for this reason that an understanding of the solid-state of a drug molecule, including the identification of all possible forms and thorough characterization of the observed forms, is of prime importance within the pharmaceutical industry. Establishing the polymorphic behaviour of a drug molecule early in development minimizes the number of unsuitable candidates progressed and reduces the risk of encountering issues later which may have major financial and time impact. Several studies can be found in literature for the characterization of mefenamic acid in different solvents and polymorphic forms obtained [2-8]. Studies are not available for the crystallisation kinetics and the determination of kinetic parameters of mefenamic acid [9], but the approach used for other chemicals can be used as a base.

Cooling batch crystallisation with seeding is used intensively to find crystallisation kinetics. These studies are interested mainly in two subjects; one is the parameter estimation of crystallisation parameters under a certain specified experimental conditions [10,11] and the other is finding the optimal cooling strategies and its effect on the final particle size distribution (PSD) of the product [12,13].

Qui and Rasmuson (1991a) used concentration and final product CSD data to determine kinetic parameters for an aqueous succinic acid system [14]. They simplified the population balance equation (PBE) to a first-order quasilinear equation system by assuming that the crystal growth rate could be expressed as a product of a function of supersaturation and a function of crystal size. They showed that by applying evaluated kinetics in process simulation allowed for a reasonably accurate prediction of the product weight mean size. 

Qui and Rasmuson (1991b) evaluated the kinetics from isothermal desupersaturation experiments by comparing polynomial fitting, initial derivatives and a nonlinear parameter estimation method [15]. In particular, this study referred to separating the volume diffusion resistance from the surface integration resistance. 

Bohlin and Rasmuson (1992) investigated the influence of growth rate dispersion on the product size distribution of batch cooling crystallisation by computer simulations [16]. The study includes three different growth rate activity distributions. For unseeded and seeded processes, the influence of the corresponding coefficient of variation was analyzed. They showed that the effect of growth rate dispersion on the crystal size distribution was significant at moderate dispersion. At high dispersion, even the actual shape of the growth rate activity distribution was important.

Doki et al. (1999) studied the effect of seeding on batch crystallisation using a natural cooling mode [17]. They showed that a critical seed concentration can be found to improve the product mean size by the ideal growth of seed crystals and reducing the secondary nucleation.

Mohameed et al. (2002) investigated three different cooling rate effects on the crystallisation kinetics parameters and on the final mean particle size of potassium chloride [18]. They used the gravimetric method to measure the concentration of potassium chloride in the liquid phase as a function of time. A non-linear optimisation method was applied to estimate the growth and nucleation rate parameters from the unseeded batch cooling crystallisation experiments. Their study suggested that unseeded batch crystallisation can produce reasonable mean particle sizes if an appropriate cooling rate is chosen.

Rohani et al. (2002) determined nucleation and growth kinetics of ammonium sulphate by using a laboratory type cooling batch seeded crystallizer, a rigorous mathematical model of the process, and a nonlinear optimisation program [19]. Supersaturation and crystal size distribution were employed as the major outputs of the model in the optimisation objective function. The sizeindependent growth kinetic was shown to lack the necessary parameters to accurately fit the experimental data.

Mohameed et al. (2003) developed a model-based optimisation strategy for the optimal cooling profile of batch crystallisation [20]. They used the numerical method of lines to solve the batch crystallizer model by discretization of population density function over the size length.

Hu et al. (2004) proposed a new optimisation-based identification approach for general batch cooling seeded crystallisation of ammonium sulphate to overcome insufficient information on the CSD based on simplified population balance models such as the moment equations [21]. The final-time CSD was directly used for identification. Consequently, size-dependent growth rates can be considered.

The determination of the two-step crystal growth rate parameters by isothermal desupersaturation experiments relies upon the adequate use of experiment data. The growth rate vs supersaturation relation is always almost linear on a log-log diagram over a reasonable range of supersaturation, regardless of the mechanism involved. Only within the slight curvature of the graph, the information available can be found to separate the volume diffusion resistance from the surface integration resistance. Although using a different polynomial to correlate the desupersaturation data actually destroys essential information and due to significant end affects, the range over which reliable growth rates may be obtained is reduced. Direct estimation of growth rate parameters by fitting the supersaturation balance equation to the supersaturation measurements circumvents these deficiencies.

Simple, one-dimensional crystallization models that neglect crystal agglomeration and breakage are common in the literature and provide insight into the inner workings of crystallization processes in terms of growth and nucleation kinetics. These models allow for the calculation of crystal size distribution outcomes and are often used in process optimization routines to determine an optimum trajectory (temperature, anti-solvent, supersaturation) that will result in a predefined crystal size distribution [22].

Liu et al. (2004) chose ciprofloxacin hydrochloride crystallisation and developed a mathematical model for the calculation of the linear crystal growth rate in solution crystallisation [23]. They used transient CSD data, the solution concentration curve, and the suspension density of a batch crystallisation process without taking into consideration the secondary processes. An advantage of this method was that crystal shape factors did not effect the linear crystal growth rate.

Hu et al. (2005) used the final product CSD, the reactor temperature profile, and the solute concentration profile for the estimation of the kinetic parameters of crystal nucleation and growth [24]. They solved the PBE directly instead of reducing it to a set of ordinary differential equations (ODEs) or by changing the PBE to simple forms under some assumptions, so it could be applied to cases with a size-dependent growth rate and the determination of optimal temperature profiles for the batch crystallizer was discussed.

Mohameed et al. (2007), presented a methodological frame-work for the optimisation of the para-xylene crystallisation process using a developed rigorous model and model parameters estimation [25]. The kinetic model parameters were estimated by reconciling the experimental data of the para-xylene concentration and total mass of para-xylene suspended crystals with model predictions using a nonlinear optimisation technique. For this purpose, the numerical method of lines (MOL) was used to solve the model’s mathematical equations. The effect of cooling strategy on the estimated parameters was studied. To find the optimum cooling profile, a new objective function was proposed and also, a nonlinear optimisation technique was applied. The new proposed objective function maximizes the yield of crystals which have a certain required average size or greater. This objective function requires the calculation of mass crystals along the size range. In the study, the homogenous wellmixed cooling batch crystallizer was used for the paraxylene crystallisation in order to investigate its advantages. They stated that, the batch crystallizer is the suitable equipment for the study of nucleation and growth rates.

In the present study, the experiments for the measurement of crystal growth kinetics were carried out by the desupersaturation curve principle based on the measurement of the solution concentration versus time in a seeded isothermal batch experiment. To estimate the growth kinetics from the desupersaturation curve obtained, the Population Balance Modelling method was used and solved by the method of moments. The simple model developed by Garside et al. (1982) was applied for obtaining the initial values [1].

1.1. Crystallisation Kinetics

Crystallisation is a two step process of nucleation and crystal growth requiring a change of free energy (Gibbs). Both nucleation and growth of crystals depend on the degree of supersaturation, but usually to different “orders”, with nucleation exhibiting the stronger dependence and higher free energy requirement. Crystals can subsequently be disrupted or can agglomerate, thereby further changing both their particle size and numbers.

Nuclei are the first formed embryos, possibly only a few nanometers in size, which subsequently grow to produce tangible crystals. Nucleation is thus the first formation of the solid phase. It occurs due to clustering and aggregation of molecules or ions in a supersaturated melt, solution or vapour to a size at which such entities become viable i.e. they will grow rather than redissolve. Crystal growth is a diffusion and integration process, modified by the effect of the solid surfaces on which it occurs. Solute molecules/ions reach the growing faces of a crystal by diffusion through the liquid phase. At the surface, they must become organized into the space lattice through an adsorbed layer. Neither the diffusion step nor interfacial step, however, will proceed unless the solution is supersaturated. The rate of crystal growth can be expressed as the rate of displacement of a given crystal surface in the direction perpendicular to the face. Different crystallographic faces of a crystal, however, usually have different linear growth rates. Variations occur in the shape of crystals when individual faces grow at different rates, the overall crystal habit being determined by the slowest growing face [26]. The crystal growth rate can also be determined by the overall growth of the crystal. In this case, the usual method to be used is to determine the change of the crystal mass per unit time [27]. The theories are classified into two groups, that is, based on either thermodynamic or kinetics properties.

1.2. Crystallisation Kinetics Correlations

Formation of new crystals may result from any one or a combination of different mechanisms (primary homogeneous or heterogeneous and secondary nucleation) or by attrition. In an industrial crystallisation practice secondary nucleation usually has a predominant influence [28]. Secondary nucleation is a complex phenomenon and is not well understood as a general theory for the prediction of nucleation rates does not exist. Several correlations based on the power law model have been found to explain most of the experimental data satisfactorily. The power law (from the Becker-Doering relationship) is given in Equation (1).

(1)

The nucleation rate constant k_b may be a function of many other variables, in particular temperature, hydrodynamics, presence of impurities and perhaps some crystal properties, and the effect of these might be incorporated in the rate equation using suitable functional forms. The power law term represents the kth moment of the crystal size distribution present in the crystallizer.

The kinetics for the secondary nucleation can be measured either by measuring the width of the metastable zone, the induction time or by counting the number of nuclei formed. One of the methods for the determination of nucleation rates is by measuring the maximum possible supercooling rate that can be obtained in a saturated solution when it is cooled at different rates [29].

The rate of nucleation can also be determined by observing the time elapsed between the creation of supersaturation and the formation of a new phase. This time interval is defined as the induction time and is a function of the solution temperature and supersaturation. The formation of a new phase can be detected in several different ways, for example, by the appearance of crystals or by changes in the properties (turbidity, refractive index) of the solution [30].

The development and operation of the industrial crystallisation processes can be made significantly easier if some data on the kinetics of crystal growth were available. This information can be incorporated in process models, can be used in process and crystallizer design, and can shed light on the observed behaviour of the system.

The crystal growth theories provide a theoretical basis for the correlation of experimental crystal growth data and the determination of kinetic parameters from the data to be used in models of industrial crystallisation processes. For engineering purposes in crystallizer design and assessment, the simple empirical power law relationship expressing the rate of mass deposition is given in Equation (2).

(2)

It is more convenient to express the overall linear growth rate instead of mass deposition as given in Equation (2) as

(3)

The constants in Equations (2) and (3) can be related to each other through the expression

(4)

In general, the overall rate constants K_g and k_g depend on temperature, crystal size, hydrodynamic situation and presence of impurities and are usually fit to the Arrhenius equation to obtain a general expression for growth rate as a function of temperature.

1.3. Determination of Crystallisation Kinetics

The kinetic parameters can be highly sensitive to small concentrations of contaminating chemicals, which can result in kinetic parameters that vary over time. Also, many crystals are sufficiently fragile that the crystals break after formation, or the crystals can agglomerate or have erosion or other surface effects that are difficult to characterize. Another significant source of uncertainty in industrial crystallizers is associated with mixing. Although crystallisation models usually assume perfect mixing, this assumption is rarely true for an Industrialscale crystallizer.

Numerous techniques have been devised to measure and analyze crystal growth and, to a far lesser extent, nucleation. The experimental techniques developed to extract crystallisation kinetic parameters illustrated in Figure 1 are based on phenomenological models.

The classical approach to establish crystallisation kinetics is to isolate the growth and nucleation processes and determine their kinetics separately by direct and or indirect methods under different hydrodynamic conditions. Direct methods to characterize crystal growth include single crystal studies; measurement of weight gain or movement of a distribution of seeds suspended in a supersaturated solution. In an indirect method, supersaturation depletion is measured.

Since all of the direct measurement techniques are time consuming and require a significant number of experiments to obtain sufficient data to obtain kinetics parameters, the indirect methods is more reasonable for estimation to both growth and nucleation kinetics. Most of the indirect methods are based on the measurement of the solution concentration versus time in a seeded isothermal batch experiment. This is often called the desupersaturation curve since the concentration and solubility can be used to calculate the supersaturation of the system versus time.

The growth rate kinetics can be estimated by a number of methods from the desupersaturation curve obtained during batch-seeded isothermal experiments. The simplest of these methods was developed by Garside et al. (1982) and involves using the first two derivatives of a second order polynomial fitting of an initial portion of a desupersaturation curve at time zero [1]. It is assumed in the method that the concentration change is due only to crystal growth with no nucleation occurring. Using this analysis, the power g and constant k_g of the crystal growth rate equation of Equation (2) can be obtained from the Equations of (5) and (6).

(5)

(6)

where ΔC₀ is the supersaturation at time zero, and are the first and second derivatives of the desupersaturation curve at time zero, respectively. F is the shape factor ratio and equals β/α. represents the average size of the seeds, A_T0 denotes the surface area of the seeds at time zero. The derivatives and are usually obtained by fitting the desupersaturation curve to the 2^nd order polynomial form of Equation (7)

(7)

So that, , and.

However simple, results obtained from this method are quite sensitive to measurement errors since the method relies on derivatives. More sophisticated techniques for the estimation of growth kinetics involve the use of the entire desupersaturation curve with parameter estimation techniques [14,15].

1.4. The Population Balance Equation (PBE) Modeling Approach

The population balance provides the mathematical framework incorporating expressions for the various crystal formation, aggregation and disruption mechanisms to predict the final particle size distribution.

The estimation technique is based on a model for seeded and batch cooling crystallisation and optimisation method which minimizes the difference between experimental data and those predicted by the model. The mathematical model comprises a population balance equation for crystal size distribution and a mass balance equation which relates the total mass of crystal size distribution and a mass balance equation which relates the total mass of crystals to the concentration change of the solution.

For a perfectly mixed batch crystallizer of constant volume, in which crystal breakage, agglomeration and crystal growth rate dispersion are assumed to be negligible, PBE for size independent growth rate is expressed as:

(8)

Subject to the boundary and initial conditions:

(9)

(10)

where n_s denotes the population density distribution of the seeds, which are added at t = 0, and G⁰ is the growth rate of nuclei.

The mass balance models the solute concentration in the continuous phase. Since the solute concentration is one of the few parameters of crystallisation that may be measured experimentally, mass balance provides a connection between the model and experiment that may be used to identify kinetic model parameters for a specific system.

The mass balance is

(11)

where C is the mass of solute per mass of solvent, r_c is the crystal density, k_v is the volume shape factor converting L³ to actual crystal volume, and h is a constant factor converting solvent mass to slurry volume.

The simplified form of the energy balance for the batch crystallizer is

(12)

where C_p is the specific heat of the slurry, ∆H_c is the heat of crystallisation, U is the overall heat transfer coefficient, V is the volume of crystallizer, A is the surface area available for heat transfer, and T_f(t) is the time dependent temperature of coolant in the jacket. Energy balance was not considered because the system is isothermal and no temperature change was observed.

The other constitutive equations required are the nucleation rate and the growth rate given as Equations (13) and (14):

(13)

(14)

with these equations the model is completely defined. The parameters k_b, b, k_g, g are inferred from experimental data.

The PBE 8 is solved using the method of moments. The i^th moment of the PBE is defined as

(15)

The first four moments are required for the general batch crystallizer model.

(16)

(17)

(18)

(19)

The initial conditions are

for i = 0, 1, 2, 3(20)

C(0) = C₀(21)

T(0) = T₀ (22)

The Equations (11), (16), (17), (18) and (19) can be solved using an ordinary differential solver. The source of stiffness in this system comes from the large difference between the kinetic parameters. The identified parameters (g, k_g, b and k_b) will be optimized using an optimisation algorithm.

2. Experimental

Mefenamic acid obtained from the Aldrich Chemical Co. with purities > 99% was used in the present study. For the estimation of the crystal growth kinetic parameters of MA recrystallized in ethyl acetate (EA), the desupersaturation curve was obtained. The solubility data required is given in Table 1.

The desupersaturation technique involves the preparation of a saturated solution of known concentration in a batch crystallizer or other vessel where the temperature can be controlled and the vessel should be equipped with a stirrer. With the stirrer on, the temperature of the solution is lowered (5˚C) to obtain the supersaturated solution.

The solution must remain clear with no crystals present at the lower temperature. Seed crystals of the solute of known mass and size or size distribution are added to the solution. Since the solution is supersaturated, these crystals will grow, causing a decline in the concentration of the solution that is measured as a function of time. This is done by taking samples and using one of a number of concentration measurement techniques. The density meter was the concentration measurement method that we used.

Figure 2 shows the apparatus to measure the desupersaturation curve online. The solution is continuously pumped from the crystallizer to a ceramic cross flow filter. Clear solution is obtained in a small side stream that returned to the crystallizer [30]. The direct measurement of the crystal weight change enables the overall growth rate to be determined as in Equation (2) at a particular temperature, crystal size and supersaturation level. This may be correlated in terms of the operating variables.

A 250 ml saturated MA-EA solution was prepared at 30˚C in a batch vessel equipped with a jacket to control the temperature and the vessel was also equipped with a stirrer. By stirring, a clear solution was obtained. After that, the temperature of the solution was lowered to 25˚C to make the solution supersaturated. The solution remained clear with no crystals present at the lower temperature and a supersaturated solution was obtained. Size distribution of the seed crystals were measured by Lasentec FBRM for calculations and added to the solution. 1% seed means that 0.0034 g of seed crystal was used. Since the solution was supersaturated, crystals grew. So, a decline in the concentration of the solution was measured as a function of time.

Samples of 2 ml were taken every minute throughout the experiment. These samples were filtered and the concentration change of the solution was measured by a densitometer. It was very difficult to measure the density because of the manual feeding of the samples into the machine. Air bubble occurrence had a negative affect on the measured density values, so the sample should be introduced without any air bubbles and intensive care was taken for achieving this result.

The density determination was based on measuring the period of oscillation of the vibrating U-shaped sample tube, which was filled with the sample liquid. The measuring principle of the instrument was based on the change of the natural frequency of a hollow oscillator when filled with different liquids or gases. For the determination of the densitometer constant, k, in Equation (23), two calibration measurements of samples of know density were necessary. Equation (23) represents the difference of densities of the two samples as follow:

(23)

For estimation of density of an unknown sample Equation (23) is arranged as Equation (24).

(24)

To determine the densitometer constant, k, measurements were done at 30˚C with air and distilled water. After a sufficient time for the proper temperature equilibrium, the period TP was read for the sample tube filled with air. Double distilled water was now introduced by means of a syringe into the lower entrance part of the sample tube. The filling of the sample tube can be observed when the illumination for the sample tube is turn on. Temperature equilibrium could be conveniently followed over several read out cycles; the TP-value becomes constant ±1, indicating complete temperature equilibrium.

TP values obtained for air and water and by using Equation (23), the k value was calculated. The calculations are given below.

During the concentration measurements, some unreasonable results were obtained. This may be caused by the failure in introduction of the sample or the evaporation of the solvent. The solvent evaporation affects the results significantly, but it was assumed that no solvent evaporation occurred. Therefore, the measurements were repeated until the results were satisfactory.

For the density meter measurements, firstly, calibretion should be accomplished. For this reason, the three best measurements used for the calibration curve are given in Table 2. The calculations are given below. Table 3 includes the data required for calculations. This table also shows the literature data for further calculation for the determination of crystallisation kinetics of mefenamic acid. Table 4 has the measured and calculated values for drawing the calibration curve. The results of densities calculated are given in Table 5. Figure 3(b) represents the calibration curve of the densitometer.

To calculate the densitometer constant, k, the following calculations were done.

After obtaining the constant, k, for the calibration curve, three sample solutions were prepared with the designated amounts given in Table 2 and by using Equation (24) the solution densities were calculated as follows.

By using the density values in Table 5, to convert the density to the concentration, the calibration curve was drawn as Figure 3(b).

At the end of the experiment, a concentration vs. time curve was drawn and labled as the desupersaturation curve. Figure 3(a) shows the desupersaturation curve of MA crystals recrystallized from EA and the data is given in Table 6.

3. Results and Discussion

3.1. Garside Solution

As mentioned previously, the simple Garside Method can be used for the rough estimations of the crystallisation growth kinetic parameters, the crystal growth rate constant and power of the mefenamic acid crystallisation process. The effect of the growth rate power on the crystallisation growth rates can also be observed. By using the calibration curve (see Figure 3(b)), the decreasing concentrations of MA-EA solution were calculated and given in Table 6 previously. The desupersaturation curve was drawn and a second order polynomial equation was fitted to the experimental data as illustrated in Figure 3(a).

For the Garside solution, the crystal growth rate power, g, was assumed as 1, 2, 2.5, 3, and 3.5 and all calculations were repeated using MathCAD. Firstly, the experimental concentration (see Table 6) values were introduced to the program. A second order line equation as a function of time was defined as Equation (25) and for the best fitting equation, an optimisation problem was solved by using MathCAD. This function was expressed as Equation (26) and this objective function was minimized.

(25)

(26)

Initial guess for parameters,

(27)

And a₀, a₁, and a₂ values were calculated by using MathCAD.

Now, the supersaturation at time zero and the first and second derivatives at time zero can be calculated:

After that the values were calculated by assuming different g values. For each assumption, the calculated values of, the growth rate constants, K_g and k_g are given in Table 7. In the following sample calculations, the growth rate constants K_g and k_g calculations for g equal to 1 can be seen.

The effect of g on linear and mass growth rates of mefenamic acid crystals are illustrated in Figures 4(a) and (b), respectively. For linear growth rate, while the g value was increasing, there was a decrease in the linear growth rate. Also during crystallisation time, the linear growth rate was decreasing because of the decrease in the concentration difference.

In the mass growth rate, in 5 min. intervals, while the g value was increasing, the mass growth rate is also increasing but after 5 min., the behaviour changed oppositely. It means that when the g value is increased, a decrease can be seen in the mass growth rate.

3.2. The Pbe Modelling Approach for Seeded and Unseeded Crystallisation

Population balance equations (PBEs) are encountered in several scientific and engineering disciplines, which model complex processes where the accurate prediction of the dispersed phase play a major role for the overall behaviour of the system. Some examples of these types of processes can be found in systems such as precipitation, crystallisation, aerosol, bubbly, droplet flows, and so on. The population balance provides the mathematical framework incorporating expressions for the various crystal formation, aggregation and disruption mechanisms to predict the final particle size distribution.

The estimation technique is based on a model for seeded and batch cooling crystallisation and optimisation method which minimizes the difference between experimental data and those predicted by the model. The PBE can be solved by different methods; in this study, the method of moments was applied. For simulations, MathCAD and MATLAB packages were used.

The mass balance equation (Equation (11)) was defined as Equation (28) and the moments of equation of PBE (Equations (16)-(19)) were introduced in MathCAD and MATLAB as Equations (29)-(32). Since the experiment was carried out in an isothermal batch vessel, in model identification the energy balance equation was not considered in the solution. Equations (33) and (34) represent the initial and boundary conditions.

(28)

(29)

(30)

(31)

(32)

In this study, the PBE modelling approach was applied for seeded and unseeded crystallisation. The zero moment of the PBE reflects the nucleation rate. For seeded crystallisation, the nucleation was not considered. The zero moment was set to a high value for representing the nucleation (e.g. for g = 1, μ₀ = 2.05 × 10⁵). The initial and boundary conditions for seeded crystallisation is,

(33)

(34)

and the initial and boundary conditions for unseeded crystallisation, as,

(35)

with the saturation concentration at 25˚C as 0.009888g_sol/ g_solv.

The objective function of the optimisation problem was defined in Equation (36):

(36)

For the MATLAB 7.0 simulations, several initial value problem solvers were used for the ODE model and these solvers are given in Table 8. For the selection of the best algorithm, a search was carried out and the effect of time step chosen and the error calculated between the model and the experimental data (Equation (36)) can be found in Table 9 Ode23s which is Rosenbrock method was selected for the further calculations that give the minimum

errors calculated. As in the Garside solution, the growth rate power, g, was assumed as 1, 2, 2.5, 3 and 3.5 and the same growth rate constants were used. For seeded crystallisation, the nucleation was eliminated and the best agreement was predicted for the zero moment of the PBE. But for unseeded crystallisation, the nucleation was considered and by changing the nucleation rate constant k_b and power b, the best prediction was searched. The calculated kinetic parameters for mefenamic acid by MathCAD and the same results achieved by MATLAB are given in Table 10 for seeded and unseeded crystallisation of mefenamic acid. The first four moments of the PBE of mefenamic acid crystallisation calculated for unseeded crystallisation can be observed in Figure 5 graphically.

The experimental data was obtained for the seeded crystallisation. From the behaviour of the zero moment curve, the time elapsed between the unseeded and seeded crystallisation was observed as 5 minutes and this nucleation time was added to the model as seen in Equation (37).

(37)

As seen in Figure 6, for g = 2.5, a good agreement between the modelling predictions of Garside and PBE and the experimental data can be observed. It can be concluded that, using the rough estimations of the growth kinetic parameters (g and k_g), all kinetic parameters can be calculated more precisely by PBE modelling.

3.3. Optimum Kinetic Parameters

The modelling data predicted by both methods and the experience gained can now serve as good initial values for the prediction of the optimal kinetic parameters. The Fminsearch algorithm was used in MATLAB 7.0 for the solution of the optimisation problem and predicted results are given in Table 11. The agreements can be seen in Figures 7(a) and (b) for seeded and unseeded crystallisation of mefenamic acid in EA respectively.

4. Conclusions

The development and operation of industrial crystallisation processes can be made significantly easier if some data on the kinetics of crystal growth are available. This information can be incorporated in process models, can be used in process and crystallizer design, and can shed light on the observed behaviour of the system. The estimation of the crystallisation kinetics parameters is very important for controlling crystallisation and obtaining desired products with suitable size, shape and polymorphic form. To determine the kinetic parameters, two methods were chosen, the desupersaturation curve technique (Garside model) and the PBE modelling approach. Because both methods require the desupersaturation curve, a concentration versus time profile during crystallisation and the experiments to obtain the concentration profile were carried out.

Mefenamic acid is crystallized in ethyl acetate solvent in an isothermal batch crystallizer with seed crystals. The density of the sample taken every minute is measured by the densitometer. For the densitometer, the device constant was firstly calculated by calibration of the machine with air and distilled water as 0.316589. By using the calibration curve obtained, the concentration values are calculated from the density measurements and a profile as a function of time is drawn. Firstly, the desupersaturation curve technique with the solution of the Garside method for the estimation of the growth kinetics can be used for quick estimation of kinetic parameters to design the process and modelling purposes, because the experiments are fast and the data is easy to obtain, but they do not supply very accurate kinetics data. Also it is very sensitive to experimental conditions and errors. In the Garside solution, g values are assumed as 1, 2, 2.5, 3 and 3.5 and values and the mass growth rate constants, K_g are calculated. These values are a set of good initial estimates for the PBE modelling approach. The second method is the PBE modelling approach for seeded and unseeded crystallisation of mefenamic acid. The parameters of the nucleation and growth rates have been successfully estimated for Mefenamic Acid using non-linear optimisation. MathCAD and MATLAB packages were used for the simulations and several algorithms, explicit methods and implicit methods, were applied for the solution of the ODE set but not a significant difference observed. Both the Garside and the PBE modelling approach were compared and it can be concluded that both methods are in good agreement with the experimental data as desired. The parameters estimated from both methods, then, serve initial values for the optimisation problem and the optimum kinetic parameters for unseeded and seeded mefenamic acid crystallization in a batch crystallizer are obtained. The effect of agglomeration and breakage on the PBE model will be researched as the next step.

5. Acknowledgements

The experimental work of this research was carried out by my master student Sevilay Gökbel at the Illinois Institute of Technology, Laboratory of Crystallisation, under the supervision of Prof. Dr. Allan S. Myerson. I am very grateful to Sevilay Gökbel and Prof. Dr. Allan S. Myerson.

REFERENCES

Nomenclature

Symbols

NOTES

^*Corresponding author.