Simulation of CaCO3 Crystal Growth in Multiphase Reaction

Calcium carbonate formation and aggregation processes have been studied from many years and there are widely described in the literature (i.e. Kitano et al., 1962, MontesHemandez et al., 2007, Reddy & Nancollas, 1976). However, the mechanism of the process, which depends on the way of reaction conducting (i.e. Dindore et al., 2005, Feng et al., 2007, Jung et al., 2005, Schlomach et al., 2006) is till now not fully understood and investigated due to increasing application of CaCO3 in commercial production of new materials, pharmaceuticals and many others. Calcium carbonate occurs in nature in three polymorphic modifications: rhombohedral calcite, orthorhombic aragonite usually with needle-like morphology and hexagonal vaterite with spherical morphology. The most needed from the practical point of view is the most stable thermodynamically calcite. One of its important applications area is connected with fabrication of functional solids where the fully controlled precipitation process must be applied. The big interest in this field is due to the fact that application of produced materials is determined by many strictly defined parameters. In recent years many researchers deal with application of organic additives (Bandyopadhyaya, 2001) as a template to produce inorganic materials and conduction of reaction in a macroor microemulsions or in sol-gel matrixes. Such methods give an opportunity to control of precipitation process or to modify product properties but unsolved problem remains purity of the obtained powder. There are also many ways of CaCO3 precipitation conducting without any additives (i.e. Cafiero et al., 2002, Sohnel & Mullin, 1982, Rigopoulos & Jones, 2003a). Although a lot of investigations described in the literature (i.e. Chakraborty & Bhatia, 1995, Chen et al., 2000, Cheng et al., 2004) there are still many questions about the full mechanism of crystals nucleation and growth of freshly precipitated particles. Generally, crystallization is a particle formation process by which molecules in solution or vapor are transformed into a solid phase of regular lattice structure, which is reflected on the external faces. Crystallization may be further described as a self-assembly molecular building process. So, crystallographic and molecular factors are thus very important in affecting the shape, purity and structure of crystals (Colfen & Antonietti, 2005, Collier &

. There are two established mechanisms of crystals growth described in literature (Jones et al., 2005, Judat & Kind, 2004, Spanos & Koutsoukos, 1998).The Ostwald ripening involves the larger crystals formation from smaller crystals which have higher solubility than larger ones (the smaller crystals act as fuel for the growth of bigger crystals).Another important growth mechanism revealed in recent years is nonclassical crystallization mechanism by aggregation, i.e. coalescence of initially stabilized nanocrystals which grow together and form one bigger particle (Judat & Kind, 2004, Myerson, 1999).There are some papers dealing with calcium carbonate formation through oriented aggregation of nanocrystals (Collier & Hounslow, 1999, Myerson, 1999, Wang et al., 2006) or through self assembled aggregation of nanometric crystallites followed by a fast recrystalization process (Judat & Kind, 2004).This way of particles formation control to fabricate ordered structures is inspired by processes observed in biological systems and is one of top topics of modern colloid and materials chemistry (Judat & Kind, 2004, Myerson, 1999, Wang et al., 2006).Each way of reaction conduction needs its own modeling.In the literature there are many different models and simulations (i.e.Bandyopadhyaya et al., 2001, Hostomsky & Jones, 1991, Malkaj et al., 2004,) done for different particular reactions.Each model describing the crystal formation (i.e.Quigley & Roger, 2008, Tobias & Klein, 1996, Wachi & Jones, 1991) has to take into account both particulate crystal characteristics and fluid-particle transport processes.The crystals formation and further solid-liquid separation of particulate crystals from solution involves suspension and sedimentation.During these processes solid matter may change phase from liquid to solid or vice versa.New particles may be generated and existing ones can be lost.Thus, both the liquid and solid phases are subject to the physical laws of change: conservation of mass and flow.The crystals may be also separated from fluids by flow through reactor.So, any model welldescribing the particular crystallization process has to take into account the conditions of reaction leading in the reactor.

Modelling and simulation
The behaviour of real crystallization processes is determined by the interaction of multiple process phenomena, which all have to be modelled to fully describe the process.Over the past decades simulation has become a standard tool for solution of these model equations.Different tools have been developed to solve many typical chemical engineering problems particularly for standard fluid phase processes.Also for more complex processes, as population balance models for crystallization processes (Randolph & Larson, 1988) and computational fluid dynamics (Ferziger & Perić, 1996) problems based on the Navier-Stokes equations, commercial simulation tools are available.However, still not every kind of crystallization process models can be solved with the available tools.Moreover, every particular crystallization process needs a specific treatment taking into account process parameters, hydrodynamic conditions, crystallizer construction, etc.A separate problem, which also has to be considered is connected with the accuracy of the simulation.The accurate calculations are time consuming and the accuracy is strongly connected with the way in which the simulation is performed.So, for each accurate simulation of a particular crystallization process it is necessary to elaborate both the appropriate physico-chemical description of the process and the proper way of simulation performance.
Simulation of CaCO 3 Crystal Growth in Multiphase Reaction 557

Conservation equations -Computational fluid dynamics
Conservation relates to accounting for flows of heat, mass or momentum (mainly fluid flow) through control volumes within vessels and pipes.This leads to the formation of conservation equations which enables to predict results of operation performed in defined equipment.In continuum mechanics the general equation for all conservation laws can be expressed in the following form (Spiegelman, 2004): where:  -any quantity (in units of stuff per unit volume), F -the flux of  in the absence of fluid transport, V -velocity of transport, H -a source or sink of .
To derive an equation for conservation of mass it is necessary to substitute  = (density -the amount of mass per unit volume), F = 0 (mass flux can be only change due to transport) and H = 0 (mass cannot be created or destroyed).After that we get the following equation, called, the continuity equation: In the case of conservation of energy (heat) in a single phase material (Spiegelman, 2004), the amount of heat per unit volume is  = c P T where c P is the specific heat at constant pressure (energy per unit mass per degree Kelvin) and T is the temperature.The heat flux consists of two components due to conduction (the heat flux is F = −k  T where k is the thermal conductivity, "-" because heat flows from hot to cold) and transport ( c P T V).Heat can be also created in a investigated volume due to viscous dissipation, radioactive decay, shear heating etc. (H -all the terms creating heat).Thus the conservation of heat equation assumes the form: (3) For constant c P and k, after introducing thermal diffusivity = k/( c P ), the equation can be rewritten in the form: Conservation of momentum (or force balance) can be derived (Spiegelman, 2004) assuming that the amount of momentum per unit volume is  = V and the forces which can change the momentum are connected with the stress that acts on the surface (F = -) and gravity (H = g, where g -terrestrial acceleration).So, the equation has the following form: (5) which can rewritten in the following way: www.intechopen.com For an isotropic incompressible fluids the above equation can be rewritten into Navier-Stokes equation: where: is the dynamic viscosity, P -fluid pressure.
In the case of crystallization a further conservation equation is required to account for particle numbers.This is the population balance.It is another transport equation which, based on particle formation (nucleation, growth, agglomeration, breakage, etc.), allows for prediction of particle size distribution, n(L, t) in the defined crystallizers.Generally the population balance equation (PBE) can be written in the following form (Jones et al., 2005) where: B 0 is nucleation rate, B and D are the "Birth" and "Death" functions for agglomeration and breakage of crystals, where v i is the "internal" velocity and v e is the "external" velocity.The "internal" velocity describes the change of particle characteristic, e.g. its size, volume or composition, and the "external" velocity, the fluid velocity, in the crystallizer.The "internal" velocity, for well-mixed systems, is approximated by the crystal growth rate (G): and usually assumes the following form: where: L is a crystal size.
For non well-mixed systems, the velocity derivatives, in addition to crystal growth, have to be included to the equation.The population balance is a partial integro-differential equation that can be normally solved by numerical methods, except for some simplified cases.Different numerical discretization schemes for solution of the population balance (Kumar & Ramkrishna, 1996, Nicmanis & Hounslow, 1998, Ramkrishna, 2000) and compute correction factors in order to preserve total mass are widely described in the literature (Hostomsky & Jones, 1991, Rigopoulos & Jones, 2003a, Wojcik & Jones, 1998, Wuklow et al., 2001).Computation Fluid Dynamics, (CFD) is the numerical analysis and solution of system involving all transport processes via computer simulation (Jones et al., 2005,).It is strongly dependent on the development of computer related technologies and on the advancement of our understanding and solving of ordinary and partial differential equations.Direct numerical solving of complex flows in real conditions requires a huge amount of computational power and is very much dependent on the physical models applied.That is why, an ideal model applied for such calculations should introduce minimum amount of complexity into the model equations, while capturing the essence of the relevant physics.
One of the most important flow phenomena is turbulence.If it is present in a certain flow it appears to be the dominant over all other flow phenomena.That is why successful modelling of turbulence greatly increases the quality of numerical simulations.Although, all analytical and semi-analytical solutions of simple flow cases were solved at the end of 1940s, there are still many open questions on modelling of turbulence and properties of turbulence it-self.Till now, no universal turbulence model exists yet.
In the case of crystallization, CFD involves the numerical solution of conservation continuity, momentum and energy equations coupled with constitutive laws of rate (kinetic) processes together with the population balance accounting the solid particles formed and destroyed during crystallization.So, the CFD model solution comprises both the flow properties and a particle size distribution what leads to the formation of conservation equations which enables to predict results of operation performed in defined equipment.
Attempts to generate a theoretical model-based description of the interaction of fluid dynamics and crystallization face the multi-scale nature of this interaction.Usually, the population balance is represented by a partial differential equation of particle size and time and the mass balance, in most cases, is expressed as ordinary differential equations.On the other hand, the growth and nucleation kinetics of particles are often based on empirical correlations.
The main problem connected with a numerical simulation is a problem of discretization of the all coordinates (Euclidean space, particle size, time).Discretization significantly affects the accuracy, the computational costs and even convergence properties of numerical algorithms.Therefore, the selection of the proper discretization grids has to be carefully considered in the context of the characteristic scales of the modelled phenomena.
Usually, for the fluid flow calculation, the Euclidean space is divided into a number of CFD grid cells with elementary volumes.The size of these cells is above the Kolmogorov turbulence scale (order of magnitude 10 -4 m) but small enough to well resolve the convective flows and energy transport within the unit (Ferziger & Perić, 1996).Such discretization is sufficient to resolve the most of the phenomena occurring in mass crystallization but needs t o b e i m p r o v e d i n t h e c a s e o f r e a c t i v e crystallization processes where micromixing phenomena play the significant role.The time coordinate of the CFD problem is also discretized using small time steps (seconds) to resolve fast fluctuations.The particle size coordinate (the population balance equation) has to be also discretized.
There are many methods available to perform this discretization (Ramkrishna, 2000, Hounslow, 1990, Hounslow et al., 1988, Hill & Ng, 1995) but in all cases, the most important in the proper evaluation of the size of CFD cell with the appropriate number of particles.If the CFD cells are too small or have too low number of particle the statistical requirements of the population balance is not fulfilled.This may result in an incorrect solution.The next problem connected with the discretization is necessity of solution of several dozens of the equations in each CFD grid cell what would certainly result in prohibitive computational cost and possibly introduce convergence problems.Therefore, some means of model reduction must be employed to allow a numerical simulation.Moreover, all these methods have been developed with a focus on the way in which the systems are mixed.Generally, CFD models can be implemented to "well-mixed" and "non well-mixed" systems.Assumption of well-mixing is commonly used for the modelling of crystallization processes, what simplifies the simulation and reduces its time.Such approach can be accepted in the case of theoretical calculations and small, laboratory scale, cristallizers.However, even in a stirred tank with impellor (Rielly & Marquis, 2001) we deal with very inhomogeneous fluid mechanical environment.The turbulence quantities and the relevant mean-flow may vary by orders of magnitude throughout the vessel, especially around the impellor.Therefore it is clear, that the 'well-mixed' assumption will lead to significant errors on the rates of growth, nucleation and agglomeration, and consequently, on the crystal size distribution.In these cases, information of the solid concentration distribution, as well as local velocities, shear rates and energy dissipation rates would be needed for the proper design of the process.Crystallization systems frequently show also high levels of supersaturation around the points where it is generated (cooling surfaces, evaporation interfaces, etc.) causing as well as suspension significant local density variations (Sha & Palosaari, 2000).Consequently, crystallization rates locally vary throughout the crystallizer even in case when no reactive crystallization occurs.Therefore, assumption of uniform conditions throughout the reactor volume can not be accepted.Moreover, many crystallization processes are directly affected by the local fluid dynamic state.One of the most important factors is the shear rate which strongly influences both the frequency and the efficiency of particle collisions (agglomeration (Hounslow et al., 2001)) as well as particle-impeller collisions (Gahn & Mersmann, 1999) which are depended on the relative velocity of the particle and local streamlines around the impeller blade.
The mixing problem increases with increasing of the scale of operation.Typically, fluid dynamics phenomena act on 'micro-scale' (CDF grid), a much smaller scale compared to crystallization phenomena which are usually considered on the 'macroscale' (unit).To solve the problem one can compute the population balance in each CFD grid, accounting for the full locality of the crystallization kinetics or use the scale or spatial resolution for the population balance.The first approach is not recommended because it can violate the statistical assumptions used for the formulation of the population balance equation and needs a long, tome consuming calculations (computational costs).The second method enables for selection of some compartments, representing a certain region in the crystallizer, which can be treated as homogenous (well-mixed) and well described by CFD.Such approach can be a compromise between one single, well mixed unit and the over-detailed system.The use of this model requires the exchange of information between the two scales of calculations: "inner" inside the compartment and "outer" between the compartments.
Taking into account these assumptions some compartmental mixing models (Wei & Garside, 1997) for modelling precipitation processes based on the engulfment theory (Baldyga & Bourne, 1984a, 1984b, 1984c) has been elaborated.Also, several mesomixing and micromixing models have been proposed to describe the influence of mixing on chemical reactions on the meso-and molecular scale (Villermaux & Falk, 1994, Baldyga et al., 1995).

Batch reactor
Batch (or semi-batch) reactor is one of the most popular reactors widely used in chemical and pharmaceutical (especially batch crystallizers) industry.Batch crystallization processes, commonly investigated, are still not well understood because the process is strongly influenced by fluid mixing, particle aggregation and particle breakage.For a batch crystallizer with nucleation, aggregation, breakage and growth occurring (Wan & Ring, 2006).the population balance equation is given by (Randolph & Larson, 1988) as: where: n(v) is the number-based population of particles in the crystallizer being a function of the particle volume v, G(v) is the volume dependent growth rate, b(v) is the volume dependent birth rate and d(v) is the volume dependent death rate.For the initial condition, n(v,t=0) = n o (v).For aggregation, the birth b a (v) and death d a (v) rate terms can be given by (Hulburt & Katz, 1964): where (v,u) is the aggregation rate constant (a measure of the frequency of collision of particles of size v with those of size u).
In the case of breakage, the birth b b (v) and death d b (v) rate terms can be given by (Prasher, 1987): where S(v) is the breakage rate that is a function of particle size v, (v,w) is the daughter distribution function defined as the probability that a fragment of a particle of size w will appear at size v.The population balance equations can be solved by the use of the standard method of moments (SMOM) and the quadrature method of moments (QMOM) (Wan & Ring, 2006).Using these methods the population balance can be simplified into a series of a few discrete moment equations (some of them as number of particles ( v m 0 ), volume of particles ( v m 1 ), etc. have physical significance) defined, for k-th volume-dependent moment, in the following way: Such calculated v m 0 and v m 1 represent the total number and total volume of particles in the system Because in the CFD code, the particle density function is described as a function of particle size x, instead of particle volume v and the population balance is written in terms of n(x) instead of n(v) the population balance, eq.[11], can be rewritten as: .
and the k-th length-dependent moment as: The both SMOM and QMOM models has been tested (Wan & Ring, 2006) using numerical cases with nucleation, growth, aggregation and breakage and the obtained results have been compared with the analytical measurements.For all cases the OMOM model gave the very good (the accuracy < 1 %) description of the particle size distribution in the batch reactor.
The particle size distribution in a batch crystallizer can be also simulated in different way.
As an example can be given a process of obtaining of calcium carbonate (Kangwook et al., 2002) when we deal with the following overall precipitation reaction: where the feeds are a solution of sodium carbonate and a solution of calcium hydroxide at certain, defined concentrations, and the main product is calcium carbonate.The main variable which is to be estimated is particle size distribution of precipitated CaCO 3 .The precipitation occurs, when the calcium ions and carbonate ions are present at supersaturated concentration levels.Supersaturation implies that the ionized species are present in the solution where the solubility of the species is exceeded.If we assume that the ionization reactions are fast compared to the precipitation i.e. the ionization reactions reach equilibrium instantaneously and that the perfect mixing in the reactor is obtained we can write the mass balance of the precipitation reactor as follows: dV F = q -q dt (19) where: C j F is the concentration of species j in the j-th feed stream, C j is the reactor concentration of species j, q j F is the feed flow rate of stream j, q F is the total feed flow rate, q is the total outlet flow rate, V is the volume of contents in the reactor, k a is the area factor, L is the characteristic particle size, G(L, t) is the growth rate of particle, n(L, t) is the particle size distribution (number of particles per volume of solvent per particle size), t is the time of reaction, j is equal to 1 for Ca(OH) 2 or 2 for Na 2 CO 3 .The mass balance equation should be solved together with the population balance equation: where: P(L.n.t) is a number of density with corresponding initial condition n(L,t 0 ) = n t0 (L) and boundary condition n(L 0 ,t) = n L0 (t) (i.e. the number of nucleated particles), where L 0 is the nucleated particle size.
Next important equation needed is an equation describing nucleation rate.Typically, nucleation and growth rates of precipitation and crystallization processes are represented by semi-empirical power laws.A proper, nucleation model has to take into account the both primary nucleation induced by supersaturation without particles and secondary nucleation related to the existing particles in the reactor.Growth rate is a function of supersaturation and particle size and can be calculated from the following equations (Eek et al., 1995): where: C s is the supersaturation of the solute, which is defined as: [(C 1 C 2 ) 0.5 −1] (expressed in terms of the normalized concentration) and an a n , b n , a t , b t , a L , b L are the parameters.
The kinetic equations have strong nonlinearity due to the power terms what combined with the mass balance equation makes the problem difficult.However, the computationally demanding part of the precipitation reactor model is the population balance equation.In general, the population balance equation can be converting into a set of ordinary differential equations.Many, various forms of the finite element method and the finite difference method can be applied for this purpose.The details on the solution techniques can be found in a (Ramkrishna, 2000) book.
The population balance equation can be simplified in the case when the right-hand side of eq. ( 20) is a linear or an independent function of the density number.Then a closed-form of the solution can be obtained using the method of characteristics (Varma & Morbidelli, 1997).
We can further simplify the model equations assuming that the aggregation and breakage are negligible and the growth rate takes a separable form of G t G L in eq. ( 22) (where: G t is the time-dependent part of the growth rate and G L is the size-dependent part).In this case we get the following equations for the particle size distribution: where t 0 is the start time of growth reaction, L b in the birth size and t b is the birth time of the L size particle at time t, which can be obtained by solving the following equations: In the case of size dependent growth, there are no general theoretical kinetics and the separable form is the exclusively used empirical form.
In order to simulate the precipitation reactor, the mass balance and the population balance equation should be solved together.They can be solved used an explicit integration method in which the algebraic equations are solved just once at the beginning of each integration step and held constant or finite element method (Kangwook et al., 2002).
The usefulness of this model for the calcium carbonate precipitation (both an explicit integration and finite element method gave almost the same results) has been checked successfully by (Kangwook et al., 2002) but is necessary to remember, that the assumption of negligible agglomeration and breakage (limits of the model) can be applied only for the reactor where the particle density is maintained on the low level (Kataki & Tsuge, 1990).
For the calcium carbonate precipitation in the batch reactor, breakage can be treated as a negligible phenomenon but the agglomeration is usually significant according to the high particle density in the reactor (Collier & Hounslow, 1999).So, if we want to avoid the aggregation and breakage phenomena in this reactor we have to operate the process in a special way, maintaining the low particle density.
Generally, the presented approach can be implemented for simple precipitation reaction and is, especially, very useful in the case of "run-to-run" or "on-line" controlling of the particle size distribution in a batch (or semi-batch) reactors (Kangwook et al., 2002).

Crystallization in tube
Every model describing crystallization in a tube has to take into account the fluid dynamics, the fluid flow through the tube and crystallization processes acting simultaneously.The simplest model describing crystallization from solution with feed concentration c 0 , in a wallcooled tube with a defined length and radius, where the supersaturation is generated by cooling of the solution by means of an energy withdrawal at the wall, can be derived making the following assumptions (Kulikov et al., 2005): the system is considered to be quasi-homogeneous -it is assumed that the flow through the tube causes very well mixing of the fluid and solid (very small crystals) phases.So, instead of writing separate transport equations for the fluid and the solid phases, a single equation for the whole suspension is formulated.This results in assuming no slip and no particle drag which also implies no segregation of the particles, -mixture properties (density , molecular viscosity , specific heat capacity c p , thermal conductivity ) are assumed to be constant, -no heat of crystallization is released, -agglomeration and particle breakage are not considered.The fluid dynamics of the homogeneous mixture can be described by the Reynoldsaveraged Navier-Stokes equations consisting of the equations for mass and momentum conservation: where:  V is the vector of Reynolds-averaged velocities, p is the static pressure, g is the gravitational acceleration, and t are viscosity and the turbulent viscosity, respectively.Using the introduced assumptions, a boundary condition was set for the flow at the tube inlet, no-slip condition was used at the wall and the standard k-model (Ferziger & Perić, 1996) has been used, as a turbulence model for a closure of the system.So, the energy balance can be expressed as: where: T is the temperature and a boundary temperature condition is specified at the walls of the tube.
The population balance equation used in this model has been taken from (Marchisio et al., 2003).It contains (Kulikov et al., 2005) the accumulation term, the particle growth term, the convective transport term, terms reflecting molecular and turbulent diffusion of particles with the molecular diffusion coefficient D m and the turbulent diffusion coefficient D t , respectively, as well as particle birth b and death d terms and can be written in the following form: where: n(L, x, t) is a particle size distribution and L is a characteristic particle size.
In this case, it is assumed a simple kinetics with the growth term G obeying McCabe's law (size-independent growth) and being first order dependent on supersaturation: where: c and c s is the solution concentration and the equilibrium concentration at saturation, respectively.k 1 is a constant.As it assumed both birth term B and death term D are set to zero: and nucleation B 0 is accounted for as a left boundary condition as follows: and expressed by a power law equation: where: is the volume fraction of solids and k 2 and k 3 are constants.The initial condition for nucleation are given by the following equations; and mass balance for the solute in the liquid phase by the following: where: D c is the solute diffusivity, cr is the density of the crystals and k v is are the shape factor of the particles.
The presented model specified in eqs ( 27)-( 36) is a multidimensional dynamic problem containing partial differential equations formulated in spatial coordinates x, one internal particle size distribution coordinate L and the time coordinate t.The locally distributed velocities, temperatures, and particle size distribution are the unknown variables which cannot be calculated analytically and have to be obtained by a numerical simulation.
As it was mentioned before the numerical simulation can be done using two approaches.
The first aims at the reduction of the complexity of the population balance discretization by selection a small number of variables characterizing the particle size distribution.It causes some loss of accuracy in the solution of the population balance, which is reformulated in terms of these variables.Transport equations are also reformulated for these variables and solved along with the CFD problem on the proper spatial grid.Usually, these variables are the moments of the distribution function i.e. the Quadrature Method of Moments (Marchisio et al., 2003).A main disadvantage of this approach is the inaccurate reconstruction of the particle size distribution when no a-priori information about its shape is available.
The second approach is based on the reduction of the spatial resolution for the population balance only.Most crystallization phenomena like growth, agglomeration, etc. do not change significantly on the resolution of the CFD grid and can be considered to act on larger scales.This allows for the representation of the population balance by collecting a set of CFD cells in an 'ideally-mixed' compartment.The population balance equations can then be solved in this compartment by a highly accurate discretization scheme.Set of such ideallymixed compartments represents different regions of the crystallizer.This approach has been well described in the literature (Kramer et al., 2000).
It is difficult to claim the superiority of one of these approaches over the other.The proper selection of the approach very much depends on the application to which it is addressed.The compartmental approach better describes the major crystallization phenomena in a cooling crystallizer with complex breakage and aggregation behaviour while the reduced population balance approach better describes a high spatial fluctuation of supersaturation, e.g., in reactive crystallization.

Bubble column reactor
Bubble column reactor is an apparatus in which simplicity of design gives rise to extraordinary complexity in the physical and chemical phenomena.That is why modeling of the precipitation process in this reactor needs an integration of reaction kinetics, population balance and hydrodynamic principles.Such successful modeling of the bubble column reactor applied for the precipitation of calcium carbonate by carbon dioxide absorption into lime has been done by (Rigopoulos. & Jones, 2003a).They used their own (Rigopoulos. & Jones, 2003b) finite element method for solving the time-dependent population balance equation with combined nucleation, growth, agglomeration, and breakage.The previous studies of gas-liquid precipitation (Rigopoulos. & Jones, 2001) which used the method of moments, took into account only a nucleation growth.However, experiment in both gas-liquid (Wachi & Jones, 1991) and liquid-liquid (Tai & Chen, 1995, Collier & Hounslow, 1999) precipitation of CaCO 3 have evidenced the presence of agglomeration and demonstrated its importance in determining of the product crystal size distribution.The time-dependent population balance equation (Rigopoulos. & Jones, 2003b) that describes the evolution of the particle size distribution in a finite, spatially uniform domain, with particle volume as the "internal" coordinate, and including nucleation, growth, and agglomeration, can be written as follows: where: n(v,t) and n in (v,t) is the population density at the reactor and at the inlet, respectively; G(v) is the volumetric growth rate; and B 0 , a , and V 0 are the nucleation rate, agglomeration kernel, and volume of the nuclei, respectively,  is a width of boundary layer.
The equation should be solved with the following initial and boundary conditions: The mass balance equation is derived from the concept of penetration theory (Astarita, 1967) where mass transfer and chemical reactions at the interface are treated simultaneously.The interface and the bulk are considered as two separate dynamic reactors which operate independently and interact at discrete time intervals.Thus, the diffusion and reaction of chemical components is described by the following equations (c i is the concentration of component i, superscripts I and B denotes variables at the interface and bulk, respectively): with initial and boundary conditions: ,0 0 In the most cases of the bubble column the precipitation phenomena at the interface can be neglected because of the very short contact time between the reagents compared to the bulk.Generally.nonideal mixing should be considered in the column but for the relatively short height of the column and intense recirculation a full mixing in the bulk can be assumed.In this case the mass balance in the bulk can be described in the following way: where f j (c i B ,n j ) is a function into which the original population balance is transformed via the finite element discretization.Initial conditions are calculated from the mixing of bulk and interface at the end of the previous contact time.The solution of the interface equations is obtained numerically with an implicit iterative scheme, while the bulk equations are calculated according to Adams method (Hindmarsh, 1983).
The reactions occurring during the process can be described in the following way: CO 2(g) + OH -= HCO 3 - ( 47) The first step (eg.( 46)) is a CO 2 absorption in water at the gas-liquid equilibrium.The equilibrium can be described by Henry's law, taking into account that we deal with an ionic system, as follows: where: H and H 0 is a Henry's constant for an ionic and nonionic system, respectively, I i is an ionic strength of component "i" and h i , h -, h + , h g is a component "i", anions, cations and gas contribution, respectively.The kinetics of carbon dioxide absorption into alkali solutions are determined by the conversion of CO 2 (aq) into HCO 3 -(eq.( 47)), which proceeds at a great, but finite rate.This reaction is followed by an instantaneous ionic reaction eq.( 48). and the precipitation reaction eq.( 49).The rate of CaCO 3 (s) production is determined by a crystallization mechanism but always volumetric crystal growth is size-dependent even when the linear growth is size-independent (McCabe's law).To obtain the rate of change for the whole crystal mass is necessary to integrate the volumetric growth function over the whole range of crystal volumes: where: CaCO3 and M CaCO3 is calcium carbonate density and molar mass, respectively.
To estimate the rate of crystal mass production, which is coupled with the population balance, it is necessary to derive a complete kinetic model of precipitation taking into account the whole information concerning the crystal formation i.e. nucleation, crystal growth as well as agglomeration and breakage.
The growth rate kinetics is usually described by the linear growth rate (the increase in particle diameter or radius) G l , by the following expression: where: k g is a kinetic constant and s is the saturation ratio (Rigopoulos. & Jones, 2003b).
Nucleation process can have a variety of mechanisms (homogeneous, heterogeneous, secondary, etc).In the bubble column it can be assumed (Rigopoulos. & Jones, 2003b) that in the beginning of the process high supersaturation levels induce primary nucleation, but later, secondary nucleation causes the rise of crystal growth.So, the overall nucleation model consists of the sum of the two models: primary and secondary.The primary nucleation which depends mainly on supersaturation is usually described by a power law: Secondary nucleation is induced by the existing crystals (Garside & Davey, 1980) and is a function of the crystal mass (M c ): where: k n1 , k n2 , k n3 , k n4 and k ns are the appropriate constants (Rigopoulos. & Jones, 2003b).Thus, the overall nucleation model can be expressed as follows: It is necessary to point out that calcium carbonate can appear, during the precipitation process, in three different polymorphs where the most prevailing polymorph appears to be calcite.That is why, a kinetic model should account for their simultaneous presence in the solution (Chakraborty & Bhatia, 1996).Usually, because of the complexity and difficulty of such calculations, the considerations are limited only to calcite.Agglomeration of crystals is a very complex and system-dependent process.Usually, it can be simplified and treated as a two-step process.The first step of agglomeration, i.e. the formation of flocculates through collisions and interparticle attraction, is similar to the phenomena occurring in colloids and aerosols.The second step is the growth of crystalline material between the clusters at so-called cementing sites (Hounslow et al., 2001).In the case of the bubble column (Rigopoulos. & Jones, 2003b) the agglomeration can be assumed to be roughly proportional to growth (Hounslow et al., 2001) and described as the second-order dependence on supersaturation, in the following way: where: k a is the agglomeration constant (Rigopoulos. & Jones, 2003b).Hydrodynamics of the gas-liquid precipitation strongly depends on the gas holdup which determines the rates of the chemical phenomena.The Eulerian-Eulerian multiphase CFD model (Rigopoulos, & Jones, 2001), where the turbulence in the liquid phase is calculated with k-model (Schwarz & Turner, 1988), can be used for its description.This model can be www.intechopen.com Modern Aspects of Bulk Crystal and Thin Film Preparation 570 successfully used for modeling large-scale equipment because it gives sufficiently accurate results with respect to averaged properties.However, it is less successful in reproducing fine details i.e. the radial phase distribution.
The above model very well (good agreement with the experiment) described the precipitation of CaCO 3 by CO 2 absorption into lime, in the bubble column.The conjunction of penetration theory and CFD predictions of the gas holdup seems to yield an adequate description of the reactor performance.Such integration of the population balance, reaction kinetics and hydrodynamic principles allowed for proper formulation of modeling approach for the gas-liquid precipitation process and the model can be used as a tool for the analysis and scale-up of industrial-class equipment.

Thin film reactor
The another approach (Kędra-Królik & Gierycz, 2010) is necessary when the precipitation goes in the thin film.It happens in the Rotating Disc Precipitation Reactor (Kędra-Królik & Gierycz, 2006, Kędra-Królik & Gierycz, 2009) used for calcium carbonate production.The reaction in liquid phase goes in contact with continuously flowing gaseous carbon dioxide in the thin film formed on the surface of the rotating disc (Kędra-Królik & Gierycz, 2006, Kędra-Królik & Gierycz, 2009).This creates a constant surface area of gas-liquid interface and the carbonation reaction of lime water involves gas, liquid and solid phase.The reactions occurring during the process are described by eqs.(46-49).
The model (Kędra-Królik & Gierycz, 2010) has taken into account not only kinetics of the multiphase reaction but also crystal growth rate.The film theory (Wachi & Jones, 1991, Danckwerts, 1970) describes the mass balance of reactants in these reactions as follows: where: t is time, c CO2 , c OH , c CO3 are the concentrations of gas reactant (CO 2(g) ), liquid reactant (OH -) and the product (CO 3 2-), respectively, G', B' are rate of nucleation and crystal growth, respectively; k is second order chemical reaction constant; D CO2 , D OH , D CO3 are diffusivity of (CO 2(g) ), (OH -) and (CO 3 2-), respectively.The component (CO 3 2-) is formed by reaction (48) and consumed by the precipitation reaction (49).It is assumed also that the concentration of (CO 3 2-) is constant across the diffusion layer.Thus the population balance of the precipitated particles is given by the following equation (Hill & Ng, 1995): where: N is a population density of particles, G is linear growth rate; L is a coordinate of particle dimension; D P is the diffusivity of particles.Substituting: N = P / L, we get: where: P(x,L i , t) is a number of density discretized in L i , L i is a particle size coordinate, L 0 is an effective nucleic size, for newly nucleated particles.
In the case of the precipitation of CaCO 3 in the thin film the small crystals are obtained due to the very high nucleation rate compared to the crystal growth rate (Kędra-Królik & Gierycz, 2010).For such very small particles, the diffusivity of the crystals (D P ) within the liquid film can be described by the Stokes-Einstein equation (Hostomsky & Jones, 1991): where: k B is the Boltzmann constant, T is temperature, is viscosity and r is radius of particle The number rate of nucleation (J n ) and linear crystal growth (G) can be expressed by the Nielsen equations (Hounslow, 1990): where: n, g -the orders of nucleation and growth, respectively; c, c * -the concentration and equilibrium saturation concentration, respectively; k n , k g -nucleation and growth rate constants, respectively.The equations can be rewritten to the following forms: where: c Ca is the concentrations of (Ca 2+ ), K sp is solubility of the product (calcium carbonate).The corresponding mass based rate equations both of nucleation and growth can be expressed by the following equations (Wachi & Jones, 1991, Danckwerts, 1970): where: is crystal density.

572
The boundary conditions for the gas-liquid interface, assuming that except for the gaseous reactant (CO 2(g) ) every component is non-volatile, are as follows: at x = 0; t > 0 → c CO2 = c CO2 0 , dc OH /dx = 0, dc CO3 /dx = 0, dP/dx = 0 (70) and for the film formed on the disk surface, assuming that newly nucleated particles have an effective nucleic size equal to L 0 , as follows: Solving the mass balance equations (eqs.58-60) and population equation (eq.62) with the boundary conditions we can calculate both the discretized density number of particles (P(t,x,L i )) and discretized diameter L i .The model describes properly the change of precipitation rate in the liquid film.The Ca(OH) 2 concentration decreases because of the very high nucleation rate.Higher supersaturation leads to smaller mean crystal size, since the nucleation rate is much more sensitive to the level of supersaturation than the growth rate.It agrees very well with the experiment (Kędra-Królik & Gierycz, 2006, Kędra-Królik & Gierycz, 2009) and is caused by the fact that the high level of supersaturation is accumulated within the liquid film due to the large diffusion resistance.The proposed model very well describes the CaCO 3 crystals formation in the rotating disc reactor and can be used and recommended for accurate calculation of the particle size and distribution obtained by gas-liquid precipitation in the reactor.However, it is necessary to remember that the model has not taken into account agglomeration of the obtained crystals and cannot be used for calculation of the aggregation process in the reactor.

Conclusion
The aim of this paper was to present the different approaches to the proper and accurate modeling and simulation of CaCO 3 formation and growth in multiphase reaction.This very complex problem has been presented for most popular, different types of reactors, i.e. batch, tube and thin film reactor as well as bubble column.
The batch (or semi-batch) precipitation process has been described by closed-form solution of population balance equation, which has not taken into account aggregation and breakage, what simplifies the simulation.However, the presented strategy is general and can be applied to batch or semi-batch processes described by more complex types of population balance equations.
In the case of tube reactor integration of simulation of crystallization and fluid dynamics was successfully applied by means of the Method of Moments.The used method allowed for reconstructing the solids fraction profiles on the fine CFD grid, while preserving the full information on particle size distribution on the coarser compartment scale.The technique is well established and has moderate computational costs.The thin film reactor has been described by the model which takes into account both kinetics of the multiphase reaction and crystals growth rate.Results of calculation agreed very well with the experiment and the model described properly the change of precipitation rate from bulk liquid to the film region and showed that the higher supersaturation leads to smaller mean crystal size, since the nucleation rate is more sensitive to the level of supersaturation than the growth rate.
The gas-liquid precipitation process in the bubble column was modeled by integration the population balance, reaction kinetics and hydrodynamic principles.The used model welldescribed the precipitation of CaCO 3 by CO 2 absorption into lime and can be recommended the analysis and scale-up of industrial-class equipment.It gave also some explanations for the experimental results.It showed that the crystal mean size increase after the pH drop is due to the disappearance of the smaller crystals by dissolution, the secondary nucleation take place because a new wave of nucleation-growth is induced by the existing crystals and crystal agglomeration starts to take place at relatively high pH and proceeds to a considerable extent because the aggregates are less frequently disrupted than in stirred tanks.Moreover, a wide review of different methods and approaches to the accurate description of crystallization processes as well as main CFD problems has been presented in this chapter.It can serve as a basic material for formulation and implementation of new, accurate models describing not only multiphase crystallization processes but also any processes taking place in different chemical reactors.
Combined population balance and kinetic models, computational fluid dynamics and mixing theory enable well prediction and scale-up of crystallization and precipitation systems but it is necessary to remember that each process (performed in the well defined reactor) needs always its own modeling.