Reactive Crystallization of Brushite under Steady State and Transient

6774

Ind. Eng. Chem. Res. 2003, 42, 6774-6785

Reactive Crystallization of Brushite under Steady State and Transient Conditions: Modeling and Experiment A. Tadayyon, S. M. Arifuzzaman, and S. Rohani* Department of Chemical and Biochemical Engineering, The University of Western Ontario, London, Ontario, Canada N6A 5B9

In recent years, there has been a significant effort among researchers to develop calcium orthophosphate based biomaterials and biomedical devices for prosthetic applications and dental implants. Although the properties of these calcium orthophosphates are quite well-known, very little information is available on the growth and nucleation kinetics of these substances. In the present study, the reactive crystallization of calcium phosphate was studied in a batch system at 25 °C. The initial calcium and phosphorus concentration in the reactor was varied from 0.01 to 0.04 mol/dm3, and the solution pH was kept under 6. Analysis by XRD, FTIR, atomic absorption and UV spectroscopy revealed that the crystals precipitated from the solution were pure brushite (dicalcium phosphate dihydrate). From the experimental determination of crystal size distribution, and with the aid of a rigorous mathematical model of the process in conjunction with a nonlinear optimization program, nucleation and two-dimensional growth kinetics of brushite crystals were determined. The existence of the two-dimensional growth rate was supported by the plate shape morphology of the growing crystals studied by scanning electron microscopy (SEM). Introduction The precipitation of calcium phosphates from aqueous solution has been studied for their industrial applications in wastewater treatment, fertilizer industry, geochemistry, water softening, food and diary industry, and most prominently for biological systems. Particularly in the study of structure and solubility of bone minerals, hyroxyapatite (HAP) and brushite have been studied extensively. It has been reported that bone is formed by the infusion of an organic matrix, composed of 25 wt % collagen and 75 wt % calcium phosphate for body support and movement.1-4 Hydroxyapatite (HAP) Ca10(OH)2(PO)8, octacalcium phosphate (OCP) Ca8H2(PO4)6‚5H2O, monetite CaHPO4, and brushite CaHPO4‚2H2O are different crystalline calcium orthophosphates that have applications in biological mineralization. Moreover, formation of two amorphous phases (ACP1 and ACP2) has also been reported.5,6 Very little information is available in the literature on the mathematical modeling and the growth and nucleation kinetics of these crystalline forms from an aqueous solution. In the present study, a mathematical model was developed to describe the nucleation and crystal growth kinetics of brushite. Brushite has particular importance in the food industry for mineral enrichment of various products and in the pharmaceutical industry as a pelletizing aid and thickening agent.7 Recently, a brushite cement was developed by Landuyt et al. for reinforcement of osteosynthesis screws.8 Development of brushite bone substitute cement has also been reported.9 This cement is biocompatible, resorbable, osteoconductive, and injectable. It readily hardens in physiological conditions. Purity and crystal size distribution are the two most important features for these applications. To * To whom correspondence should be addressed. Tel.: (519) 661-4116. Fax: (519) 661-3498. E-mail: [email protected].

produce crystals with high purity and definite crystal size distribution, the kinetic study of brushite crystallization is very important. Precipitation of orthophosphates from aqueous solution is very complicated, because of the possible occurrence of many different phases. Different solid phases incorporate different ions or complexes. Concentrations of these ions change considerably with the change of calcium and phosphate concentrations and the solution pH. Boistelle and Lopez-Valero10 reported that below pH 6.5, brushite is the predominant phase, whereas above this pH the formation of ACP and octacalcium phosphate is more pronounced. On the other hand, Wu and Nancollas11 claimed that although brushite and octacalcium phosphate form readily from solution, they are thermodynamically metastable with respect to HAP and serve as a precursor phase. The order of the precipitation from a supersaturated solution is governed not only by the thermodynamic solubility product but also by kinetic factors. So, the solid-phase first precipitated undergoes changes in solution toward phases of higher stability.12 In the present study, solution pH was maintained below 6 during the experiments, and the initial molar concentration of calcium and phosphorus ranged from 0.01 to 0.04 mol/dm3. The precipitated crystals were characterized by X-ray Diffraction, Fourier Transform Infrared Spectroscopic (FT-IR) analysis, and atomic absorption (AA) and UV spectroscopy for elemental calcium and phosphorus. The XRD pattern and FT-IR spectra showed a very good match with the published information for dicalcium phosphate dihydrate. The Ca/P atomic ratio found in the precipitated crystals also matched brushite. Therefore, it was concluded that the precipitated phase was brushite. The nucleation and growth kinetics of brushite were estimated by comparing the experimental particle size distribution with the theoretical dynamic model. Three

10.1021/ie030354+ CCC: $25.00 © 2003 American Chemical Society Published on Web 11/26/2003

Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6775

Figure 1. Schematic diagram of experimental setup: 1, calcium and sodium (solution A) chloride solution flux; 2, phosphorous solution (solution B) flux; 3, NaOH solution flux; 4, turbidity sensor; 5 and 6, peristaltic pumps; 7, nitrogen tank; 8, jacketed reactor; 9, motor; 10, pH electrode; 11, thermocouple; 12, pH meter; 13, field point; 14, turbidimeter; 15, computer, V1, V2, and V3 valves.

cases were considered in modeling the process. In the first case the solutions containing calcium and phosphate ions were mixed together at a particular ionic strength without changing the solution pH. A solidphase did not form in this stage. Experimental evidence supported the validity of this assumption. Solutions containing calcium and phosphate ions were kept for 48 h, without manipulation of pH. No crystals were formed during this period. The mathematical model was able to estimate the equilibrium concentration of different species and ions in the solution. In the second case, the model was extended to estimate the concentration of different species and ions at equilibrium condition, in the presence of solid precipitation from the solution. In the final case, the dynamic reactive crystallization was considered with the aid of population balance to estimate the nucleation and growth kinetics of brushite. Experimental Section Materials and Method. In the present investigation a 3-L stirred tank crystallizer made of Pyrex glass was used. Mixing in the reactor was assured by running a propeller type impeller at 535 rpm. Nitrogen gas was bubbled into the solution to prevent the dissolution of any CO2 that might result in carbonate apatite formation. Reactive crystallization was carried out by mixing equal volumes (700 mL) of sodium phosphate (NaH2PO4‚ H2O, EM Science, NJ) and calcium chloride (CaCl2, Fisher Chemicals, NJ). The schematic diagram of the experimental setup is presented in Figure 1. All experiments were performed at room temperature. Initial molar concentration of calcium and phosphorus was varied from 0.01 to 0.04 mol/dm3 in the reactor. The feed solutions were prepared in two separate containers under nitrogen atmosphere in demineralized water. Precisely weighed sodium chloride NaCl (EM Science, NJ) was added with CaCl2 solution to maintain a fixed ionic strength of 0.15. Accumet combination pH electrode (Cole Parmer 55500-18) along with Omega DP 24 pH meter were employed to measure solution pH in the reactor. A type J thermocouple was used to measure the reactor temperature. Experiments were started by adding calcium and phosphorus solution by a Masterflex

Figure 2. Typical variation of pH and % transmittance with time.

peristaltic pump (L/S 7520-40, Cole Parmer Instrument Co.) to the reactor. A representative sample stream from the crystallizer was continuously circulated through the sensing zone of a turbidimeter sensor, where the transmittance was measured. The suspension was finally recirculated to the crystallizer. The sensor is composed of a simple circulation loop that includes a glass tube and an infrared emitter-detector assembly. The outputs of the turbidimeter, pH meter, and the thermocouple were sent to a computer via Field Point and Lab View (National Instruments, TX) data acquisition system. The sampling interval was 30 s. Experimental Results. A series of experiments with different calcium to phosphorus molar ratios was conducted. The initial molar concentrations used during different experiments are presented in Table 1. Once the reactants were added to the reactor, the pH of the solution reached equilibrium very quickly. The pH was then increased gradually by adding 0.1 N NaOH at a flow rate of 6.4 mL/min by a peristaltic pump. The nucleation point was assumed to be when the transmittance dropped by 1%. Because delay in recognizing the nucleation point, practically meant pumping more NaOH solution to the reactor, we had to utilize unfiltered data for detecting nucleation point. It was needed to select transmittance drop threshold in a way to ensure that the measured drop is really due to the nucleation rather than the noise. After several careful preliminary experiments, it was observed that the variations of the raw signal were close to but always less than 1%. Based on these, it was decided to select transmittance drop by 1% as the onset of the nucleation point. At this time the base addition pump was turned off. The pH was at its maximum at the nucleation point and decreased gradually due to the formation of H+ ions and then reached equilibrium. The pH and transmittance profiles of a typical experiment are presented in Figure 2. Point ‘A’ indicates the onset of base addition. Point ‘B’ represents the nucleation point. At this point the sodium hydroxide pump was turned off. Point ‘C’ represents the apparent equilibrium when pH in the reactor almost reached a steady-state value. It is notable that the percent transmittance is 99% when the pH is at its maximum value supporting the assumption that the nucleation point is achieved when transmittance drops

6776 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 Table 1. Molar Concentration Used during Different Batch Experiment molar concn of CaCl2 (mol/L) molar concn of NaH2PO4‚H2O (mol/L)

molar concn of CaCl2 (mol/L) molar concn of NaH2PO4‚H2O (mol/L)

molar concn of CaCl2 (mol/L) molar concn of NaH2PO4‚H2O (mol/L)

batch 1

batch 2

batch 3

batch 4

0.01 0.01

0.02 0.02

0.03 0.03

0.04 0.04

batch 5

batch 6

batch 7

batch 8

0.01 0.005

0.02 0.01

0.03 0.015

0.04 0.02

batch 9

batch 10

batch 11

batch 12

0.01 0.02

0.02 0.04

0.03 0.06

0.04 0.08

Figure 3. X-ray diffraction of sample.

by 1%. A sample was withdrawn from the reactor 7 min after the nucleation point was reached for particle size analysis using Mastersizer Hydro MU2000 (Malvern Instruments, UK). Another sample was drawn at the same time for magma density measurement. The particle size and magma density thus determined were used for estimating the nucleation and growth kinetics of brushite. The experiments were ended when the pH in the reactor reached an apparent steady-state condition. At the end of each experiment, the slurry was taken out of the reactor and filtered. The crystals collected on filter paper were dried in a vacuum drier. The crystal morphology was studied by scanning electron microscopy (see Figure 5). The phase identification was achieved by XRD (see Figure 3), FT-IR spectroscopy (see Figure 4) and elemental Ca and P analysis by atomic absorption spectroscopy and UV spectroscopy. Figure 3 demonstrates that the peaks corresponding to the crystals formed in our experiments match almost perfectly the standard peaks of pure brushite. Figure 4 is the FT-IR spectrum of the major functional groups of the product samples obtained in our experiments. The

Figure 4. FTIR spectra of CaHPO4‚2H2O recorded at room temperature.

fingerprints match closely those of brushite reported in FTIR standards.18 Further confirmation that the crystals obtained were indeed brushite was achieved by

Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6777 + Ca2+ + H2PO4 T CaH2PO4

K7 )

γCaH2PO+4 [CaH2PO+ 4] γCa2+ γH2PO-4 [Ca2+][H2PO4]

(7)

H2O T OH- + H+ K8 ) KW ) γOH- γH+ [OH-][H+] (8) where [ ] shows the molar concentration and γ represents the activity coefficient of various species. The latter is calculated from the Davis equation14

I1/2 log γ( ) -A|Z+Z-| - 0.3I 1 + I1/2

where Z is the species charge number and I is the ionic strength. The coefficient A depends on temperature and is calculated from

Figure 5. SEM picture of leaflet shape of brushite crystals.

Ca/P atomic ratio equal to 1. Calcium analysis was done in a Spectra AA55 atomic absorption spectrometer (VARIAN, U.S.A.), and the phosphorus analysis was done as a phosphomolybdate by absorbance measurement in a Cary 100 UV-vis spectrophotometer (VARIAN, U.S.A.). Theory Case 1: Modeling Aqueous Calcium Phosphate System without Precipitation at Equilibrium. When the two reactants are added together without adjusting pH to the nucleation point, the resulting solution reaches equilibrium with no precipitation. The reactions are proton transfer or diffusion based and exhibit rapid kinetics; therefore, they reach their equilibrium state in a short time.13 The reactions and their associated reaction equilibrium constants are + H3PO4(aq) T H2PO4 + H + γH2PO-4 γH+ [H2PO4 ][H ] K1 ) (1) [H3PO4]aq 2+ H2PO4 T HPO4 + H

K2 )

HPO24 T PO34 +

PO34 +

2+

Ca

T

+

H

+ γHPO2γH+ [HPO24 ][H ] 4

γH2PO-4 [H2PO4] K3 )

(2)

+ γH+ [PO3γPO34 ][H ] 4

γHPO2[HPO24 ] 4

(3)

CaPO4 K4 )

γCaPO-4 [CaPO4] γCa2+ γPO3[Ca2+][PO34 ] 4

Ca2+ + OH- T CaOH+ K5 )

γCaOH+ [CaOH+] γCa2+ γOH- [Ca2+][OH-]

(4)

[CaHPO4]aq 2+

γCa2+ γHPO2[Ca 4

][HPO24 ]

A ) 0.486 - 6.07 × 10-4T + 6.43 × 10-6T2 (10) K1 to K8 are equilibrium constants and can be found in ref 15:

K1 ) 10-2.148, K2 ) 10-7.199, K3 ) 10-12.35, K4 ) 106.4, K5 ) 101.3, K6 ) 102.74, K7 ) 101.4, K8 ) 10-13.997 In eq 9, Z and I are known, and therefore the activity coefficient, γ, can be calculated. The material balance equations are + [Ca]T ) [Ca2+] + [CaPO4 ] + [CaOH ] +

[CaHPO4] + [CaH2PO+ 4 ] (11) 23[P]T ) [H3PO4]aq + [H2PO4 ] + [HPO4 ] + [PO4 ] + + [CaPO4 ] + [CaHPO4] + [CaH2PO4 ] (12)

The system of equations has 13 variables and 10 equations (eqs 1-8, 11, and 12). [Ca]T and [P]T are the known initial concentrations of calcium and phosphorus, and pH of the solution is fixed or measured by the experimenter, therefore, is a known value. This sets the system with 10 equations and 10 unknowns. The algebraic system of equations is highly nonlinear; a mathematical solution for this system will be discussed in a subsequent section. With the set of equations described above, pH remains as a degree of freedom; this allows the experimenter to change the pH by the addition of a base and enter the value in the model as an input. If one desires to theoretically find the equilibrium pH of the solution when no base is added, then a charge balance needs to be added to the above set of equations: + 2[Ca2+] + [CaOH+] + [CaH2PO+ 4 ] + [H ] )

(5)

23[H2PO4 ] + 2[HPO4 ] + 3[PO4 ] + [OH ] +

[CaPO4 ] (13)

Ca2+ + HPO24 T CaHPO4 K6 )

(9)

(6)

In this study, a charge balance was not included in the equations in order to have the freedom to manipulate the pH.

6778 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003

Case 2: Modeling Calcium Phosphate System with Precipitation at Equilibrium. As mentioned earlier, the solution resulting from the addition of two reactants does not lead to the precipitation of crystals if the pH is not adjusted to the nucleation point. When the pH is increased, depending on the initial reactant concentrations, at a certain pH, the solution becomes supersaturated with respect to calcium phosphate and triggers nucleation. If the pH of the solution is not controlled, the reactions shift toward the direction of generating more H+, which in turn reduces the pH of the solution. The process is fast in early stages after the nucleation but gradually reaches an apparent equilibrium. Depending on the pH of the process, various phases of the product may be formed. For pH less than 6.0, brushite is thermodynamically the most stable form, while at pH greater than 6.0, hydroxyapatite is the final phase of the product. For ionic reactions, it is known that as the product species concentration exceeds the solubility product, Ksp, the system becomes metastable with respect to the compound and the substance precipitates. In the case of brushite, this is expressed as16 -6.622 [Ca2+] [HPO2) 4 ] g Ksp (Ksp ) 10

(14)

To model the precipitation process at equilibrium, eqs 11 and 12 need to be modified to reflect the precipitation of brushite + [Ca]T ) [Ca2+] + [CaPO4 ] + [CaOH ] +

[CaHPO4] + [CaH2PO+ 4 ] + [CaHPO4‚2H2O]s (15) 23[P]T ) [H3PO4]aq + [H2PO4 ] + [HPO4 ] + [PO4 ] + + [CaPO4 ] + [CaHPO4] + [CaH2PO4 ] + [CaHPO4‚2H2O]s (16)

where the subscript “s” indicates the precipitated solid phase. Eqs 1-8, 14 (taking ‘equality’ corresponding to the onset of nucleation), 15, and 16 are the major constituents of the model for equilibrium precipitation. Case 3: Dynamic Modeling of Reactive Crystallization (Precipitation). A model of reactive crystallization is composed of two major parts: the reactions that contain the calculation of various species concentrations in the solution and the crystallization that includes the computation of the dynamic size distribution and magma density. A dynamic model requires the dynamic mass balance equations meaning that reaction equilibrium constants (Keq) should be replaced with the forward and backward reaction constants (kf and kb). For example, reaction 3 will be k3f

+ 8 PO3HPO24 79 4 + H k 3b

(17)

This adds seven pairs of reaction constants (for eqs 1-7) to the system of equations. Estimation of the reaction kinetic parameters (k1f, k1b, k2f, k2b, ....k7f, k7b) requires many experiments and solving difficult optimization problems. Fortunately, it is known that these reactions are fast and equilibrium is established instantaneously; hence the dynamics of the reactions may be ignored. On the other hand, the nucleation and the growth of calcium phosphate are relatively slow with the corresponding time constants in the order of min-

utes; therefore, the dynamics of the precipitation process need to be taken into account. This means that once the pH of the solution exceeds the nucleation point, new crystals are generated. Formation of these tiny crystals and growth of available crystals require slow integration of Ca2+ and HPO24 with each other and with the existing solid phase. This process momentarily pushes the system out of equilibrium; however, the solution quickly readjusts itself to a new equilibrium state. In the new equilibrium state, if the solution is still supersaturated with respect to calcium phosphate, solids mass is again deposited. This process is continued until the level of supersaturation becomes negligible. It is noted that in practice, mass deposition and the solution readjustment to establish the new equilibrium state are simultaneous, meaning that the concentrations of various species at any time are close to their equilibrium at the solution pH. Assumptions used in the development of the dynamic reactive crystallization are as follows: 1. The reactions are fast, the dynamics of the reactions in the aqueous solution are ignored, and equilibrium equations are used for the determination of various species in the reactions. However crystallization mechanisms including nucleation and growth are dynamic. 2. Only brushite is produced. This is accomplished by keeping the pH of the reaction below 6.0. Figures 3 (XRD) and 4 (FTIR) and elemental calcium and phosphorus analysis performed by AA and UV spectroscopy showed that the crystals are pure brushite. 3. Crystals are born at size zero. This is a common assumption in many crystallization modeling studies. 4. Crystals are plate shape. Scanning Electron Micrographs (SEM) of the samples taken at various times during batch runs in this study confirm this assumption (see Figure 5). This indicates that the horizontal (radial) and vertical (axial) growth rates are different. A shape factor is used to convert the nongeometric crystal plates to cylindrical plates (cylinders with radius r and a thickness x). At any time the population is n(r,x,t). The diameter of the cylinder is taken as the arithmetic mean of the largest (L1) and smallest (L2) dimensions of the plate like crystals, i.e., 2r ) (L1 + L2)/2. For any cylinder, there is a sphere with a diameter d that possesses the same volume as the cylinder (Vcylinder(r,x) ) πr2x ) Vsphere ) (1/6)πd3). For these cylindrical crystals, growth rate in the radial direction is much larger than that in the axial one. 5. Radial and axial growth rates are size-independent. 6. Only primary nucleation occurs. This is supported by the fact that the level of supersaturation in the reactive crystallization is extremely high compared to that in the cooling crystallizers; therefore, the solution is highly unstable and leads to primary nucleation. 7. No agglomeration or breakage is assumed. In this study, the stirrer is adjusted at a low speed to minimize the breakage of crystals. Mass Balance. The dynamic mass balance for H3PO4 is

d[H3PO4] ) -k1fγH3PO4 [H3PO4] + dt + k1bγH2PO-4 γH+ [H2PO4 ][H ] (18) Based on the assumption of fast reactions, the accumulation term d[H3PO4]/dt is set to zero. This leads to an equilibrium equation shown in eq 1. The mass

Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6779 2+ balance for [HPO24 ] and [Ca ] involves the integration of the dissolved species with the solid particles

d[HPO24 ] ) dt + γH+ [HPO2[H2PO(-k2bγHPO24 ][H ] + k2fγH2PO4 ]) + 4 4 (-k6fγHPO24

2+ γCa2+[HPO24 ][Ca ]

+

(-k3fγHPO24

(

J ) K exp

)

λ (ln(β))2

+

+ k3bγPO3γH+ [γPO3][H ]) 4 4

MW(HPO24 )

dM(CaHPO4‚2H2O) (19) dt MW(CaHPO4‚2H2O) where MW shows the molecular weight of the ionic species. In eq 19, the chemical equilibrium assumption results in + γH+ [HPO2[H2PO-k2bγHPO24 ][H ] + k2fγH2PO4] ) 0 4 4 2+ γCa2+[HPO2-k6fγHPO24 ][Ca ] + 4

k6bγCaHPO4,aq[CaHPO4]aq ) 0 [HPO2γH+ [γPO3][H+] ) 0 -k3fγHPO24 ] + k3bγPO34 4 4 (20)

(25)

The growth rate of crystals for precipitation processes may be expressed as19

Gi ) kgi(β - 1)ni

k6bγCaHPO4,aq[CaHPO4]aq) + [HPO24 ]

Assuming that the reaction takes place at a constant temperature, eq 23 may be rewritten as

(26)

where kgi and Gi are the growth rate constant and growth rate in a specific direction. The term (β - 1) in eq 26 ensures mass deposition due to growth continues until the level of supersatuartion approaches 1. K, λ, kgi, and ni are model parameters and estimated by the optimization algorithm. Population Balance. Population balance is the main equation determining the product size distribution. For batch experiments it is expressed as

∂n(r,x,t) ∂(Grn(r,x,t)) ∂(Gxn(r,x,t)) + + ) 0 (27) ∂t ∂r ∂x As discussed earlier, the population balance is twodimensional in space. This is because the growth rates are different in the radial and axial directions. Equation 27 is a hyperbolic partial differential equation, and its solution requires an appropriate numerical scheme. The boundary conditions for eq 27 are

Therefore

d[HPO2MW(HPO24 ] 4 ) )× dt MW(CaHPO4‚2H2O) dM(CaHPO4‚2H2O) (21) dt The term dM(CaHPO4‚2H2O)/dt represents the rate of change in the solids mass due to crystallization; since it is a slow process, the change in [HPO42-] is gradual and therefore dynamic. Similarly mass balance on [Ca2+] results in

MW(Ca2+) d[Ca2+] )× dt MW(CaHPO4‚2H2O) dM(CaHPO4‚2H2O) (22) dt The mass balance on the [Ca]T and [P]T (eqs 15 and 16) are still valid for the dynamic precipitation model. Crystallization kinetics includes nucleation and growth of crystals. The model for primary nucleation is4

J ) K exp

(

)

-fv2∂3 (KBT)3 (ln(β))2

(23)

where J is the nucleation rate (no./kg solvent/s), f is a shape factor, KB is the Boltzmann constant, v is the volume of a molecule in the crystal structure, ∂ is the interfacial free energy, and β is the supersaturation defined as

β)

[Ca2+][HPO2γCa2+γHPO24 ] 4 Ksp

(24)

n(0,0,t) )

J G|r ) 0, x ) 0

n(0,x,t) ) 0 n(r,0,t) ) 0

(28)

If the process is initiated in the absence of seeds, the initial condition is

n(r,x,0) ) 0

(29)

Solution of the Models. Case 1. In this case, the reactants are added together, the reactions proceed and reach equilibrium state. However as pH is below the nucleation point, no crystal is precipitated. Thirteen 2+ variables are [H3PO4]aq, [H2PO4 ], [H ], [HPO4 ], 32+ + [PO4 ], [Ca ], [CaPO4 ], [OH ], [CaOH ], [CaHPO4]aq, [CaH2PO+ 4 ], [Ca]T, and [P]T. Out of these, [Ca]T and [P]T are the initial concentrations of calcium and phosphorus and are known. [H+] is fixed by the addition of NaOH and is measured by a pH-meter. The pH usually represents the independent variable (x-axis) in the distribution diagram and is introduced to the model as a known value. This fixes the system with 10 equations (defined by 8 reaction equilibrium equations, and Ca and P mass balances) and 10 unknowns. The Davis equation for the calculation of the activity coefficients of various species does not add new unknowns to the system. By setting the temperature of the reaction and the solution ionic strength, A and γ in eqs 9 and 10 can be easily calculated. The system of algebraic equations presents eight highly nonlinear and two linear equations. Two solution algorithms are discussed: the first one follows a special procedure leading to an analytical solution, while the second solution uses the symbolic toolbox in the Matlab (Mathwork, Massachusetts, U.S.A.) software. The first algorithm is briefly


described here, more information can be found in ref 16. The algorithm starts with defining the ratio of the concentrations of various species to the total concentration of phosphorus:

R0 ) [H3PO4]aq/[P]T R1 ) [H2PO4 ]/[P]T R2 ) [HPO24 ]/[P]T R3 ) [PO34 ]/[P]T R4 ) [CaPO4 ] / [P]T R5 ) [CaHPO4]/[P]T R6 ) [CaH2PO+ 4 ]/[P]T

(30)

From eq 12 one can write

R0 + R1 + R2 + R3 + R4 + R5 + R6 ) 1 The reciprocal of R0 is

[P]T [H2PO[HPO24] 4 ] 1 ) )1+ + + R0 [H3PO4]aq [H3PO4]aq [H3PO4]aq [PO34 ] [H3PO4]aq

+

[CaPO4] [H3PO4]aq

+

[CaHPO4]aq [H3PO4]aq

+

[CaH2PO4+] [H3PO4]aq

(31)

Various terms on the right-hand side of eq 31 can be calculated using eqs 1-8. For example, rearranging eq 1 gives

[H2PO4] [H3PO4]aq

)

K1

1 γ [H ] H2PO-4 γH+ +

(32)

Equations 1 and 2 can be used to calculate [HPO24 ]/[H3PO4]aq:

[HPO24 ] [H3PO4]aq

)

K1K2

γ2

[H+]2 (γH+)2γHPO24

(33)

The same approach may be used to calculate all terms on the right-hand side of eq 31. Substituting these values in eq 31 gives the value of R0. Knowing the initial concentration of phosphorus, [P], the concentration of [H3PO4]aq can be calculated from

[H3PO4]aq ) [P]TR0

symbolic solutions become too long to be used in the computer program. In such situations, it is suggested to reduce the size of the system of equations by removing an equal number of equations and unknowns from the system. Executing the Matlab program calculates the unknowns of the system in terms of the known variables. This procedure usually leads to a solution, which in turn is a system of algebraic equations but can be easily solved with a simple numerical approach. Case 2. In this case brushite is precipitated; its molar solids concentration adds a new unknown to the system. On the other hand, the precipitation equation describing the mathematical relationship between the molar concentrations of the two ionic species of brushite at equilibrium adds a new equation. Therefore the system is composed of 11 equations and 11 unknowns and can be solved by the second approach explained in the previous section. Case 3. This case involves both reaction and crystallization phenomena, of which the latter is dynamic. The precipitation condition (eq 14) is not used because it represents the precipitation reaction under steady-state conditions. The nucleation rate (eq 23) defines the number of new crystals generated in the system per unit time per kg of the solvent. The growth rate (eq 26) sets the mass deposition of the dissolved solutes on the surface of the existing crystals, and the population balance fixes the population density of crystals at various time steps. The integration of Ca2+ and HPO24 with the solid phase is the rate-limiting step; the slow changes in the ion concentration shifts fast reactions from one side to the other side. As mentioned earlier, the population balance is twodimensional in space and one-dimensional in the time domain. Several numerical schemes for the solution of population balance exist. In this work a two-dimensional Lax-Wendroff numerical scheme is used. The derivation of the equation is presented in Appendix B. The computational size and time intervals in the simulation programs are ∆x ) ∆r ) 0.01 µm and ∆t ) 0.1 s. Estimation of Crystallization Kinetics. Out of the three cases described above only case 3 predicts size distribution and therefore uses crystal nucleation and growth correlations in its structure. The estimation of the crystallization kinetics of calcium phosphate has not been studied extensively in the literature. In this work an attempt is made to address this issue. An optimization algorithm was used to estimate the crystallization kinetics given in eqs 25 and 26, using the experimental particle size distribution. The objective function in the optimization program is to minimize the difference between the measured and the simulated percent number of crystals

(34)

The same procedure may be used to calculate the concentrations of other species. This procedure is timeconsuming, tedious, prone to calculation errors, and therefore is not recommended for the solution of the model. The second algorithm employs Matlab for seeking symbolic analytical solutions for the nonlinear algebraic equations. In this work, we used “solve” function in the symbolic toolbox for the solution of the nonlinear system of equations. The procedure is easy to learn and convenient to use. A simplified computer program for case 1 is included in Appendix A. Occasionally the

Ne

min Φ(θ) ) θ

Nc

wq[ypq,m - ypq(θ)]2 ∑ ∑ p)1 q)1

(35)

where θ denotes the kinetic parameter vector and ypq,m and ypq(θ) represent the measured and the predicted population densities, respectively. Ne is the number of experiments and Nc is the number of channels of Malvarn Mastersizer. The weight factor, wq, is assigned for each channel of the Mastersizer to recognize the relative importance of various size ranges. To optimize the crystallization kinetic parameters in eqs 25 and 26, the function fminsearch in the optimization Tool Box of Matlab was used. This function employs


Figure 6. Logarithm of the equilibrium concentration of various species vs pH ([P]T ) 0.001, [Ca]T ) [0.001]).


the Nelder-Mead Simplex approach17 for the minimization of the objective function. A simplex is a geometrical figure formed by R + 1 points in an R-dimensional space. In two-dimensional space, a simplex is a triangle; in three-dimensional space, it is a pyramid. The algorithm is based on evaluating the objective function at the vertices of a simplex, finding the worst vertex of the simplex, forming its symmetrical image through the center of the opposite face and generating a new simplex based on the new vertex. This shrinks the volume of the simplex if the newly found point is better than the older point. This step is repeated until the volume of the simplex is smaller than some user-defined preset value. Results and Discussion Case 1. In this case, the pH of the solution is lower than the nucleation point, and the reactions reach an equilibrium state without generating crystals; meaning that the nucleation condition (eq 14) is not realized. Figure 6 shows the logarithm of equilibrium [Ca2+], aqueous [CaHPO4], [HPO24 ], and logarithm of the formation product of calcium phosphate vs pH. The initial concentrations of calcium and phosphorus in the reactant solutions are low ([P]T ) 0.001 mol/kg solvent and [Ca]T ) 0.001 mol/kg solvent). The experimental observations showed that varying pH up to 9 does not lead to the nucleation of crystals. The formation of brushite occurs over the pH range from 3 to 6. In this range, the logarithm of the formation product varies between -7.6 to -10.5; these are less than log(Ksp) ) -6.22, and therefore nucleation is not expected. The fact that crystals are not formed between pH 3 and 9 indicates that the system is not supersaturated with respect to any form of calcium phosphate solids including, among others, brushite, tricalcium phosphate, amorphous, and hydroxyapatite. Figure 7 shows the 23logarithm of [H2PO4 ], [HPO4 ], [H3PO4]aq, and [PO4 ]. As it can be seen, with increasing pH the concentration 3of [H3PO4]aq and [H2PO4 ] decrease and that of [PO4 ] increases. In Figure 8, the variations of logarithm of + [CaHPO+ 4 ], [CaPO4 ], [OH ], and [CaOH ] with pH are sketched. From Figures 6-8, it is noticed that Ca2+ and H2PO4 possess the dominant ion concentrations in pH range of 3 to 6.



Case 2. Figures 9 and 10 show the equilibrium concentrations of various species in the reactor when the [Ca]T ) 0.02 and [P]T ) 0.03 molar, respectively. The figures are prepared for pH less than 6, because brushite



Figure 12. Equilibrium concentration of brushite vs pH ([P]T ) 0.02, [Ca]T ) variable; shown on the diagram).

Figure 11. Equilibrium concentration of brushite vs pH ([P]T ) 0.02, [Ca]T ) 0.03]).

Figure 13. Differential solids concentration vs pH.

crystals are generated in this pH range. As it can be seen, except for the calcium ion, the relation between logarithms of the equilibrium concentrations of various species with pH is linear. Figure 11 displays the simulated molar concentration of brushite vs pH. As it can be noticed, for pH less than 4.2, the nucleation condition is not met, and therefore no crystal is generated. Above this pH, the solids concentration (or magma density) increases with pH in a nonlinear manner and almost reaches a plateau at pH ) 6. Figure 12 sketches the equilibrium concentration of brushite vs pH, for variable [P]T concentrations between 0.02 and 0.06 molar, and [Ca]T ) 0.01 molar. At a fixed value of pH, the solids concentration increases with increasing phosphorus concentration; however, the amount of the increase is different for various pH. This is better shown in Figure 13, where the amount of crystals generated in different runs is compared to that of the first run where [P]T ) 0.02. The amount of the increase observed in the generated solids increases with pH and passes through a maximum at pH ) 4.8 and then decreases up to pH ) 6. Case 3. In this case, in addition to the mass balance, the population balance is also included in the system of equations. This enables the model to predict both the

dynamic crystal size distribution. The population balance requires nucleation and growth information to predict the outputs. The dynamic number percent distributions from 12 experiments (listed in Table 1) were provided to the optimization algorithm for the estimation of kinetic parameters. In each experiment, samples were withdrawn twice, once seven minutes after nucleation had begun and another after reaching apparent equilibrium. Both sets of data, each including population density for different classes of crystal size, were utilized to estimate the kinetics parameters. The optimizer was set to estimate the parameters based on the minimization of the error between the measured and the modeled number distributions. Malvern Mastersizer 2000 reported number percent of the particles vs size, assuming that particles were spherical, giving the distribution as a function of volumetric spherical equivalent diameter. As mentioned earlier, in this study the population density, n(r, x), was two-dimensional in space with crystals being disklike or cylindrical. To make the size distribution reported by the instrument and the model comparable, the model output was converted to a new distribution based on the spherical equivalent diameter. The optimization algorithm gave the following kinetic correlations:


More detailed analysis on the size distribution reported by the model indicates that some crystals may grow as large as 70 µm, although the spherical equivalent size is limited to 9 µm. Figure 15 shows the dynamic magma density and the supersaturation level. It is also noticed that after around 960 s, the suppuration is not fully depleted, showing that the process will continue for a longer time and the magma density in the crystallizer will increase at a lower rate. It is seen that the supersaturation ratio reduces from 78 immediately after the nucleation to almost 10 at the time of size sampling at t ) 960 s. Conclusion Figure 14. Comparison between experimental and simulated size crystal size distributions at t ) 960 s, [P]T ) 0.08, [Ca]T ) 0.08.

Figure 15. Dynamic magma density and supersaturation ratio, [P]T ) 0.08, [Ca]T ) 0.08.

(

)

-1.01 J ) 1 × 1028 exp (ln(β))2

Gx ) 1 × 10-11(β - 1)1.01 Gr ) 7 × 10-13(β - 1)0.98

(36)

Figure 14 compares the experimental and the predicted crystal size distributions with the above kinetics correlations. The predictive capability of the model was tested for the remaining 6 experiments whose size distribution were not used in the parameter estimation. Figure 14 represents the predictive capability of the model in the worst case, and it is not obtained by fitting. As it can be seen on the size axis, the equivalent spherical diameter ranges between 1 and 9 µm. The majority of crystals are born in the early stages of the process where the level of supersaturation ratio is high. These crystals grow as reactions proceed and leave little supersaturation at the end of the reaction, resulting in relatively insignificant generation of the new crystal at the end of the experiment. Based on the same graph, the number of crystals with an equivalent spherical diameter of 5.5 µm shows a maximum on the diagram. It is noted that based on eq 36, the radial growth rate of crystals is much higher than that of the axial, resulting in the generation of plate or disklike crystals.

Reactive crystallization of calcium phosphate was studied in a batch reactor at room temperature and a solution pH below 6.0. The percent transmittance and temperature were continuously monitored and recorded during the course of the reaction. Crystal size distribution of the sample crystals precipitated during the dynamic crystallization period was measured by a particle size analyzer. A mathematical model was developed to estimate the crystallization kinetics using a nonlinear optimization algorithm and the experimental CSD data from 12 experiments. The crystallization kinetics developed in this fashion were used to predict the final CSD and showed good agreement. Three of the six estimated parameters of nucleation rate and twodimensional growth kinetics, namely λ, n1 and n2 are close to unity and can be set to 1. Based on XRD, solidstate FT-IR, atomic absorption and UV spectroscopic study, it was concluded that the precipitated phase was brushite. Crystal morphology was studied by SEM. Twodimensional leaflet shape crystals were observed in this study, with a very large aspect ratio in radial and axial direction. Appendix A: Matlab Program % This simplified Matlab m.file solves a system of nonlinear algebraic equations for % modeling equilibrium reactions when two reactant (NaH2PO4 and CaCl2) solutions are % added together clear all % Define symbolic variables syms H3PO4 H2PO4_ H HPO4_2 PO4_3 Ca2 ... CaPO4 OH_ CaH2PO4 P_tot CaOH CaHPO4 ... k1 k2 k3 k4 k5 k6 k7 k8 % Write all equations and unknowns and solve the system of algebraic equations S ) solve (′k1 ) H2PO4_ * H/H3PO4′, ′k2 ) HPO4_2 * H/H2PO4_′, ... ′k3 ) PO4_3 * H/HPO4_2′, ′k4 ) CaPO4_/PO4_3 /Ca2′, ... ′k5 ) CaOH/OH_/Ca2′, ′k6 ) CaHPO4/Ca2 /HPO4_2′, ...′k7 ) CaH2PO4/Ca2/ H2PO4_′, ′P_tot ) H3PO4 + H2PO4_ ...+ HPO4_2 + PO4_3+CaPO4_+CaH2PO4+CaHPO4′,...′k8)OH_*H′, ′H3PO4, H2PO4_, HPO4_2, PO4_3, CaPO4_, ... OH_, CaH2PO4, CaOH, CaHPO4 ′) % Print results (concentration of species) ′ H3PO4 ) ′, S.H3PO4, ′ H2PO4_ ) ′, S.H2PO4_, ′ HPO4_2 ) ′, ... S.HPO4_2, ′ PO4_3 ) ′, S.PO4_3, ′ CaPO4_ ) ′, S.CaPO4_, ... ′ OH_ ) ′, S.OH_, ′CaH2PO4) ′, S.CaH2PO4, ′ CaOH ) ′, ... S.CaOH, ′ CaHPO4 ) ′, S.CaHPO4


Appendix B: Two-Dimensional Lax-Wendroff Method Population balance:

∂n ∂n ∂n + g1 + g2 ) 0 ∂t ∂x ∂r

(27)

We start with the Taylor series expansion for n(x,r,t+∆t) 2

n(x,r,t+∆t) ) n(x,r,t) +

∂n ∂2n (∆t) ∆t + 2 + O(∆t)3 ∂t 2 ∂t (B-1)

() (

)

l l ∂2n ∂ ni,j+1 - ni,j-1 ∂ ∂n ) ) ∂r∂x ∂x ∂r ∂x 2∆r l l ∂ni,j-1 ∂ni,j+1 ∂x ∂x ) 2∆r

l l l l ni+1,j+1 ni+1,j-1 - ni-1,j+1 - ni-1,j-1 2∆x 2∆x ) 2∆r l l l l ni+1,j+1 - ni-1,j+1 - ni+1,j-1 + ni-1,j-1 ) )E 4∆x∆r

2

l+1 l ) ni,j + ni,j

∂2n (∆t) ∂n ∆t + 2 + O(∆t)3 ∂t 2 ∂t

Substitute these in eq B-4:

(∆t)2 + 2 (∆t)2 + O(∆t)3 g1g2E(∆t)2 + g22D 2

l+1 l ) ni,j - g1A∆t - g2B∆t + g12C ni,j

Rearranging eq 27 gives

∂n ∂n ∂n ) -g1 - g2 ∂t ∂x ∂r

(B-2)

Taking the derivative from both sides with respect to time leads to

and rearrange to obtain the final form of the equation: l+1 l ni,j ) ni,j - ∆t(g1A + g2B) +

∂ ∂n ∂ ∂n ∂ 2n ) -g1 - g2 2 ∂t ∂x ∂t ∂r ∂t

( ) ( ) ∂ ∂n ∂ ∂n ) -g ( ) - g ( ) ∂x ∂t ∂r ∂t ∂ ∂n ∂n ∂ ∂n ∂n ) -g -g -g -g )-g -g ) ∂x ( ∂x ∂r ∂r ( ∂x ∂r 1

2

1

2

) g1

1

∂2n

∂x

2

+ 2g1g2 2

2

1

2

∂2n ∂2 n + g22 2 ∂x∂r ∂r

(B-3)

Substitute eqs B-2 and B-3 into eq B-1:

(

l+1 l ni,j ) ni,j + -g1

(

g12

)

∂n ∂n - g2 ∆t + ∂x ∂r

)

2 2 ∂2n ∂2n 2∂ n (∆t) + g + O(∆t)3 + 2g g 1 2 2 2 2 ∂x∂r 2 ∂x ∂r 2

∂n ∂n ∂2n (∆t) l+1 l ni,j ) ni,j - g1 ∆t - g2 ∆t + g12 2 + ∂x ∂r ∂x 2 2

∂2n ∂2n (∆t) (∆t)2 + g22 2 + O(∆t)3 (B-4) g1g2 ∂x∂r 2 ∂r Calculate various derivatives: l

l

∂n ni+1,j - ni-1,j ) )A ∂x 2∆x l l ∂n ni,j+1 - ni-j-1 ) )B ∂r 2∆r l

l

l

l

l

l

∂2n ni+1,j - 2ni,j + ni-1,j ) )C ∂x2 (∆x)2 ∂2n ni,j+1 - 2ni,j + ni,j-1 ) )D ∂r2 (∆r)2

(

)

g22D g12C (∆t) + g1g2E + + O(∆t)3 (B-5) 2 2 2

Notation A ) a coefficient in eq 9 defined by eq 10 (-) f ) a shape factor in eq 23 (-) g ) growth rate exponent (-) G ) growth rate (m/s) I ) ionic strength (mol/kg solvent) J ) nucleation rate (no./kg solvent/s) K ) nucleation coefficient (no./kg solvent/s) K1-K8 ) reaction equilibrium constant (variable units) kb ) backward reaction constant (variable units) KB ) Boltzmann constant ()1.38 × 10-16 erg/mol/K) kf ) forward reaction constant (variable units) kg ) growth rate coefficient (m/s) Ksp ) solubility product (for brushite in this study: mol2/ kg solvent2) M ) magma density (kg solids/kg solvent) m ) concentration (mol/kg solvent) MW ) ionic weight (kg/kg mol of the ion) n ) population density (no./kg solvent/(µm)2) Nc ) number of size channels in Malvern Mastersizer 2000 Ne ) number of experiments R ) number of dimensions in the definition of a simplex (-) T ) temperature (K) t ) time (s) w ) weight factor in eq 35 y ) variable in optimization eq, eq 35 Z ) charge number (-) Greek Letters γ ) activity coefficient (-) σ ) interfacial free energy (erg/cm2) Φ ) objective function β ) supersaturation ratio (-) θ ) vector of optimization parameters ν ) volume of a molecule in the crystal structure (Å3) R0-R6 ) ratio of the concentration of various species to the concentration of phosphorus (-) λ ) a constant in eq 25

Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6785 Subscripts aq ) aqueous i ) radial or axial direction m ) measured variable p ) index for the number of experiments q ) index for the number of channels r ) radial s ) solid T ) total or initial x ) axial

Literature Cited (1) Betts, F.; Blumenthal, N. C.; Posner, A. S. Bone Mineralization. J. Cryst. Growth 1981, 53, 63-73. (2) Posner, A. S.; Betts, F. Synthetic Amorphous Calcium Phosphate and Its Relation to Bone Mineral Structure. Acc. Chem. Res. 1975, 8, 273-281. (3) Levinskas, G. J.; Nerman, W. F. The solubility of bone mineral. I. Solubility studies of synthetic hydroxyapatite. J. Phys. Chem. 1955, 59, 164-168. (4) Mathers, N. J.; Czernuszka, J. T. Growth of hydroxyapatite on type 1 collagen. J. Mater. Sci. Lett. 1991, 10, 992-993. (5) Christoffersen, J.; Christoffersen, M. R.; Kibalczyc, W.; Andersen, F. A. A Contribution to the Understanding of the Formation of Calcium Phosphate. J. Cryst. Growth 1989, 94, 767777. (6) Christoffersen, J.; Christoffersen, M. R.; Kibalczyc, W. Apparent Solubilities of Two Amorphous Calcium Phosphates and Octacalcium Phosphates in the Temperature Range 30-42 °C. J. Cryst. Growth 1990, 106, 349-354. (7) Ullmann F. Ullman’s Encyclopedia of Industrial Chemistry, 5th ed.; Weinheim, Federal Republic of Germany, 1985; pp 498500. (8) Landuyt, P. V.; Peter, B.; Beluze, L.; Lemaitre, J. Reinforcement of Osteosynthesis Screws With Brushite Cement. Bone 1999, 25(2), 95S-98S.

(9) Charriere, E.; Terrazzoni, S.; Pittet, C.; Mordasini, P.; Dutoit, M.; Lemaitre, J.; Zysset, P. Mechanical Characterization of Brushite and Hydroxyapatite Cements. Biomaterials 2001, 22, 2937-2945. (10) Boistell, R.; Lopez-Valero, I. Growth Units and Nucleation: the Case of Calcium Phosphate. J. Cryst. Growth 1990, 102, 609-617. (11) Wu, W.; Nancollas, G. H. The Dissolution and Growth of Sparingly Soluble Inorganic Salts : A Kinetic and Surface Energy Approach. Pure Appl. Chem. 1998, 70(10), 1867-1872. (12) Abbona, F.; Madsen, H. E. L.; Boistelle, R. The Initial Phases of Calcium and Magnesium Phosphates Precipitated from Solutions of High to Medium Concentrations. J. Cryst. Growth 1986, 74, 581-590. (13) Estrin, J. Precipitation Processes. In Handbook of Industrial Crystallization; Myerson, A., Ed.; Butterworth-Heinemann: Boston, 1993; pp 131-149. (14) Buttler, J. N. Activity Coefficient and pH: In Ionic Equilibrium: Solubility and pH Calculations; Wiley: New York, 1998; pp 41-49. (15) Smith, R. M.; Martell, A. E. Critical Stability constants: Inorganic Complexes; Plenium Press: New York, 1976; Vol. 4. (16) Buttler, J. N. Polytropic Acids and Bases. In Ionic Equilibrium: Solubility and pH Calculations; Wiley: New York, 1998; pp 157-179. (17) Nelder, J. A.; Mead, R. A. Simplex Method for Function Minimization. Comput. J. 1965, 7, 308-313. (18) FTIR reinvestigation of the spectra of synthetic brushite and its partially deuterated analogues. Trpkovska, M.; Soptrajanov, B.; Malkov, P. J Mol. Struct. 1999, 480-481, 661-666. (19) Serna, H. C.; Ma, D.; Braatz, R. Optimal seeding in Batch crystallization. Can. J. Chem. Eng. 1999, June 77, p590-596.

Received for review April 25, 2003 Revised manuscript received October 3, 2003 Accepted October 13, 2003 IE030354+

Reactive Crystallization of Brushite under Steady State and Transient

Recommend Documents