Numerical investigation on freeze-drying of aqueous material frozen with pre-built pores☆

Wei Wang *,Dapeng Hu Yanqiu Pan ,Guohua Chen *

1.Introduction

Freeze-drying or lyophilisation is a dehydration method that can significantly reduce product damage caused by thermal drying.It plays a unique role in the processing of delicate and heat-sensitive materials of high values such as food,pharmaceuticals and biological products.How ever,freeze-drying is the most sophisticated and expensive process of all drying techniques,both in capital investment and in operational expenses[1].Therefore,reduction of freeze-drying time in order to lower energy consumption and raise productivity has been a world wide challenge during the past decades[2].

There are a few ways of process enhancement that can minimize drying time on a consistent basis while maintaining acceptable quality of freeze-dried products.The obvious one is to conduct freeze-drying,from the view point of energy saving,at the eutectic temperature for a liquid phase that crystallizes during freezing or the glass transition temperature for a liquid phase that remains amorphous[3].The alternative way is to combine freeze-drying with microwave heating,which has demonstrated its potential to enhance the internal heat transfer due to the interaction between the microwave and lossy media[4,5].Excessive microwave energy input,how ever,could cause frozen materials to melt and collapse,which would ruin the overall freeze-drying process.

Although modern industry has been benefited from the existing know ledge of physical and chemical processes during freeze-drying,fundamental formulations in some areas are still incomplete.Microwave freeze-drying with the assistance of dielectric materials has successfully settled the heat transfer problem.The mass transfer problem still remains unsolved.It is necessary to re-examine the freeze-drying process.There are usually three stages in freeze-drying:freezing stage,primary drying stage,and secondary drying stage.Freezing is the first stage of afreeze-drying process,and the performance of the overall process is closely related to this stage.The shape of pores,pore size distribution and pore connectivity of porous dried layer formed by ice sublimation depend largely on the size of ice crystals during the freezing stage[3,6].This dependence is of extreme importance because the mass transfer rate is significantly affected by the porous structure of the dried layer[7,8].If ice crystals are small and discontinuous,the mass transfer rate of vapour in the dried layer would be limited.However,if large ice crystals are formed and homogeneous dispersion in the frozen solution can be realized,the mass transfer rate could behigh and materials could be dried more quickly[9,10].The prevailing transfer resistance of freeze-drying is from water vapour migration in the dried region,depending on the size of ice crystals formed during freezing[11,12].Also,mass transfer parameters during freeze-drying are strongly dependent on the textural and morphological parameters of the ice phase[13].

Unlike fruits and vegetables with naturally formed porous structures that may be partially filled with water,liquid material to be freeze-dried is usually frozen into the solid material without initial voids.The latter freeze-drying operation is know n to be lengthy and expensive.It was therefore proposed that the liquid material is first frozen into the porous material with certain porosity and then freeze-dried[14].Generally,such a frozen material is partially filled with ice crystals like icecream.It would be beneficial to promoting the freeze-drying rate because of the increased area of heat and mass transfer initially.Experimental studies have show ed that the freeze-drying process can be significantly enhanced using the initially porous frozen material[15,16].The aims of the present work are to perform a numerical investigation using proper adsorption-desorption relationships to reproduce and predict the kinetics process under the experimental drying conditions,and to examine the heat and mass transfer phenomena during the freeze-drying of the initially porous material.

2.Mathematical Model

In the present simulation,a two-dimensional mathematical model of the cylindrical coordinate system was derived based on the understanding of heat and mass transfer mechanism developed previously[17].The main assumptions made in the model formulation have been reported elsewhere[2].In order to match the experimental tests,the sample geometry during simulation exactly follow ed the experimentally used one as show n in Fig.1.

In the supporting pad,

In the material region,

Fig.1.Sample geometry.

w here

Initial temperature and ice saturation are uniform.

At the central line of the sample,there are no mass and heat fluxes.

On the inner surface of the supporting pad,there is no mass flux.

On the sample surface,heat transfer is governed by radiation and drying chamber has a constant pressure.

It has to be pointed out that heat and mass flux balance equations at sublimation interface were not involved in the present model as a boundary condition.A boundary is known to be a border that encloses an area or a space.The important feature of a boundary in the sense of transport phenomena is,in the authors' understanding,that there is an abrupt change in physical properties around areas adjacent to the boundary.In the model ling of freeze-drying,an artificially moving speed of sublimation interface was usually introduced as an additional boundary condition.It is reasonable for the freeze-drying of the initially solid frozen material due to the absence of initial pores.How ever,it is difficult to justify this introduction for the initially porous frozen material based on physical considerations.The term ?R/?t is an independent term,which assumes that the moving speed of interface at any location contributes to a special heat and mass transfer that is different from other areas for the uniform,homogeneous and isotropic porous matrix.It is the additional term that results in an additional variable,i.e.,the interface position,R,in a mathematical model and causes additional complexity of model solving.Based on such understanding,some researchers presented their simulation results with a sharp change in ice saturation occurring at the interface although the materials to be dried were initially unsaturated[4,17,18],which seems to be physically unrealistic.More importantly,the sublimation interface should be inherently formed during drying,rather than inserted as a constraint for the initially unsaturated frozen material.Recent studies have con firmed this judgement[19].Current simulation results further verify that the sublimation interface can be indeed formed naturally as seen subsequently.

3.Material and Properties

During freeze-drying under arelatively low pressure,Knudsen diffusion was considered to be the controlling step for vapour diffusion in porous medium[20].

where d is the average pore diameter of porous dried product.

In a typical freezing process of aqueous solution,ice crystals have a mean size of 45-50μm[21].It has also been determined experimentally that the ice crystal size was around 50 μm[13].Therefore,50 μm of average pore diameter was used here for the initially solid frozen sample,and the value for the initially porous one can be obtained as a function of the initial saturation[19].

Based on the assumption of rigid porous medium,intrinsic porosity or the porosity of completely dried product,ε,is constant.How ever,instant aneous porosity,ε(1-S),varies with the saturation in the material region that is partially filled with ice crystals.The permeability thus increases as the moisture content decreases due to the development of additional void space.According to Kozeny-Carman equation,a permeability expression was derived based on the assumption that the dried cake texture is represented by a bundle of capillary tubes,with d being an equivalent diameter[19]as below:

It is recognized that the simple thermodynamic equilibrium relationship should be replaced by an adsorption-desorption equilibrium in the drying of hygroscopic porous materials with bound moisture[17].Unfortunately,such are lation ship has been seldom developed theoretically or experimentally due to the complexities of porous media.Therefore,an empirical relation has to be used.It is know n that the ratio of saturated vapour pressures between the adsorption equilibrium and the thermodynamic equilibrium ranges from 0 to 1[5].In other words,such a ratio should be 0 as the ending point for the absolutely dried product(S=0),and 1 as the starting point with voids filled fully with pure ice(S=1).The present simulation will test three kinds of the equilibrium relations in an attempt to examine which one can best represent the experimental results.

The first one is a simple power-law relation reported previously[5].

where P0is the saturated vapour pressure of pure ice in the thermodynamic equilibrium,which is the maximum value in the adsorption equilibrium.

The second one is a polynomial formula recommended by Redhead for multilayer adsorption[22]as

The third one is in an exponential form similar to Kelvin's equation[23].It was also constructed to satisfy the two end values of the adsorption equilibrium.

It has to be mentioned that α,β orγ in each equilibrium relation is an empirical parameter.It was determined by fitting the measured drying kinetics.Furthermore,the adsorption-desorption equilibrium was assumed always to be valid throughout the drying process.The emissivity,e,was regarded as the instantaneous surface area occupancy of solid substrate,and the angle factor,F,was considered to be unity[24].

In the present simulation,mannitol was used as the primary solute in aqueous solution and de-ionized water as the solvent.The solution employed in the formation of samples had the initial moisture content,X0of 4.48 kg·kg-1on the dry basis corresponding to 81.77 w t%on the w et basis,which are all the same in the experiments.Other physical properties can be found in Table 1.

Table 1 Physical properties used in simulation

4.Numerical Procedure

The governing Eqs.(1)-(3),together with the initial and boundary conditions(4)-(10),were solved numerically using the finite difference method.The control-volume method with the fully implicit scheme was used for discretization.The cell-centred scheme was adopted,i.e.,grid nodes were placed at the geometric centres of control volumes.The equations were discretized into a general form of linear algebraic equations for the cylindrical sample,as

w here Φ denotes the generalized independent variable,T or S.Since equation coef ficients,aP,aE,aW,aN,aSand b,strongly depend on T and S,iteration is necessary.

Because the coefficients of discretized equations typically construct a penta-diagonal matrix,the tri-diagonal matrix algorithm(TDMA)is no longer applicable.Therefore,the alternating direction implicit(ADI)method was used here.Detailed description on the ADI method can be found elsewhere[30].In fact,this method could be regarded as the alternative use of the TDMA,i.e.,a two-dimensional problem is converted into two one-dimensional implicit problems in series,each of which is solved with the TDMA.

An under-relaxation was employed to avoid the divergence in the iterative solving of such strongly nonlinear equations.Numerical errors were tested with grid numbers of 50×50,60×60 and 70×70,respectively.60×60 was finally selected because less than 1%difference was found between the predicted drying times of 60×60 and 70×70.It has to be pointed out that simulations were carried out by closely following the material geometry and experimental operating conditions.The convergence criterion of iterations is as below:

5.Results and Discussion

5.1.Comparisons between simulated and measured drying curves

In order to compare model predictions with the experimentally measured data,the initially solid(S0=1.00)and porous(S0=0.28)frozen samples were employed with the same mass and the same initial moisture content.The latter case was chosen as a representative because it has been determined to be the good operating condition from experiments[15].Sample mass was all 1.8 g with the same diameter of 14.8 mm but different heights.Detailed operating conditions are given in Table 2.Fig.2 illustrates comparisons between measured drying curves and simulated ones using different kinds of adsorption equilibrium relation.

Table 2 Tested operating conditions

Fig.2.Comparisons between experimental and predicted drying curves with different adsorption equilibrium relations.

First of all,the numerical simulation follows the trend of experimental measurement.When using the simple power-law adsorption equilibrium with the parameter α being fitted as 0.5,comparisons between the two drying curves show that the model overestimated the drying rate to a large extent in the early stage of drying.When using the Redhead type polynomial relation with the parameter β being fitted as 0.12,the model underestimated the drying rate for most of the drying process until the very end w hen the moisture content was very low.The discrepancy must be derived from the simple form of the two relations because the vapour pressure is only a function of moisture saturation.

There sults above indicate that the adsorption-desorption equilibrium relationship is a process function that involves moisture saturation and also temperature.Kelvin type equation obviously fits the experimental data the best with the parameter γ being constant at 5000.The much larger value of γ than those of α or β is because temperature effect has been taken into account in Eq.(15),with γ at the same order of magnitude as RvT.It needs to be pointed out that 5000 is a fitted value with its physical meaning yet determined.When γ is larger than 5000,the model would overestimate the drying rate,and γ less than 5000 would underestimate the drying rate.As stated previously,there exists a level of moisture content,above which the drying process is governed by the ice-vapour thermodynamic equilibrium and below which the process is controlled by adsorption equilibrium for hygroscopic materials[17].This level of moisture content has been identified as the “initial”bound moisture content[3].It depends,however,on the properties of porous material and operating conditions,namely,local temperature and pressure.Only a few properties such as intrinsic porosity,ε,and moisture saturation,S,were defined to describe a moist porous medium during drying.There is a lack of hygroscopic factor,i.e.,the affinity effect between solid matrix and moisture.Theoretically,ageneralized equilibrium relationship should be constructed to satisfy the ice-vapour thermodynamic equilibrium under high moisture content and the adsorption equilibrium under low moisture content simultaneously.More importantly,the first-order partial derivative of vapour mass concentration with respect to temperature,T,or saturation,S,is naturally continuous.Therefore,the parameter present in such an equilibrium relationship would include the hygroscopic factor.The excellent agreement achieved with the exponential form of adsorption equilibrium relation indicates that this type of equation with γ=5000 is closer to the thermodynamic equilibrium than the other tw o ones in the earlier stage of drying w hen the material moisture content is high.

It should be pointed out that the parameter,α,β or γ used in the respective equilibrium relation was kept constant for the initially solid(S0=1.00)and porous(S0=0.28)frozen materials over the w hole drying period.

5.2.Pro files of temperature and saturation within a sample

Since the exponential form of equilibrium relation gives good fitting of the experimental data,the temperature and saturation distributions are further analysed for this case under the same operating conditions aslisted in Table 2.Fig.3 show s the temperature pro files for the initially porous frozen sample(S0=0.28)during drying.The case for the initially solid one is not discussed here because it is similar to that analysed previously[5].

Initially,temperature was uniformly distributed as show n in Fig.3(a).At the beginning of drying,Fig.3(b),the temperature in the material region quickly dropped to a very low value.This is the consequence of the dynamic balance between mass and heat transfer.Instantaneously intensive sublimation would absorb a large amount of heat to balance the latent heat of sublimation at a given initial temperature driving force,leading to the highest drying rate.With the advancement of drying process,the temperature pro files,Figs.3(c-e),look similar except for the gradually increased material temperature volumetrically.At the late stage of drying,Fig.3(f-i),the temperatures on both the upper surface and the side surface increased quickly.This implies that those areas progressively went into the secondary drying stage and a relatively high temperature was required to maintain continuous internal desorption.It is important to note that the temperature close to the inner surface of the supporting pad also increased rapidly at the same time.Due to the existence of the largely pre-built pores,this part of moisture would bere moved much earlier for the porous frozen material than the solid one.This suggests that besides ambient radiation,a certain amount of heat needed for drying came from heat conduction through the pad.The pad temperature was always higher than that of the material region because of heat conduction within the pad without internal phase change.The low est temperatures were located at the cylindrical sample core,which indicates that a bulk “ice region”still existed within the sample during drying for the initially porous frozen material.The highest temperatures were predicted at the end of drying as show n in Fig.3(j).

Fig.3.Temperature pro files within a porous sample at S0=0.28 under P amb=22 Pa and T am b=30°C.

Fig.3(continued).

Fig.4 show s the ice saturation pro files within the porous frozen sample(S0=0.28)during drying.The naturally formed sublimation interface is seen in Fig.4(b-e)to recede gradually inside the sample from the upper surface and the side surface as well as the inner surface of the pad as the drying process proceeds.It is impressive to find in Fig.4(f-i)that the saturation in the area close to the inner surface of the pad decreased almost to the same extent as the other two surfaces.This phenomenon has been experimentally observed due to the initially prebuilt void space and pierced solid matrix[16].A small quantity of moisture still remained after the “interface”moving away due to the adsorption-desorption equilibrium,leaving a “tail”behind.It is notable that the saturation in the bulk “ice region”looks unchanged around 0.28,and a highland platform of saturation seems to exist always for most of the drying process.Gradually decreased saturation in this area is not observed as estimated in the previousanalysis[16].In fact,carefully examining a set of arbitrarily picked-up simulated data among Fig.4(b-h)found that there was still alittle bit decrease in saturation either along the radial direction or along the axial direction(refer to Fig.1).There was the highest value slightly above 0.28 located at the centreline of the cylinder close to the pad on the platform,and the lowest value slightly below 0.28 oriented to the upper edge of the sample surface.This implies that the platform is still an active area being subjected to dynamic sublimation de sublimation during drying.Since the adsorption equilibrium relation and the perfect gas law regulate the entire drying process,there are two dependent intensive variables among temperature,vapour concentration,pressure and ice saturation.More or less of the local saturation is closely related to the local temperature and pressure.It is known that the variable separation method was used in constructing such adsorption-desorption equilibrium relationships as stated above.With the gradually increased temperatures on the platform as seen in Fig.3(b-h),the saturated pressures,P(T,S)in the adsorption equilibrium and P0(T)in the thermodynamic equilibrium have to increase.This should attribute to a dynamic balance resultant of mass and heat transfer,leading to a similar ratio within the platform.It is expected that higher temperature and low er pressure would result in the gradually decreased saturations in this area,and vice versa.It is observed that the drying rate in the radial direction was faster than that in the axial direction throughout the entire drying process.Different from thermal drying of w et materials where the pressure difference is one of the main driving forces of liquid movement w hen liquid phase is continuous,freeze-drying has no pressure-driven flow of liquid moisture because of the solidified moisture(ice).The magnitude of ice saturation depends on the local temperature and saturation.Although saturation in the post-interface moving area was maintained at a low level,it did not yet reach the final value required for the end of drying.Therefore,additional heat should be supplied to increase the sample temperature in order to achieve the forthcoming desorption.The drying process ceased as the average residual moisture content met the product requirement and ice saturation was uniformly re-distributed,Fig.4(j).

Fig.4.Ice saturation pro files within a porous sample at S0=0.28 under P amb=22 Pa and T am b=30°C.

Fig.4(continued).

It is know n that freeze-drying of conventionally saturated frozen material is characterized by movement of sublimation interface to remove the free moisture during the primary stage,and desorption of the bound moisture during the secondary stage[3,31].For the initially unsaturated frozen material,the present simulation has demonstrated that the sublimation-de sublimation could occur throughout the entire material domain.The pre-built porous structure,pierced solid matrix and the larger internal surface area for this type of material are the reasons for the enhanced freeze drying performance with the former has much more beneficial effect on the free moisture removal by sublimation,and the latter two effects promote the bound moisture removal by desorption as show n experimentally[16].

5.3.Effect of initial porosity

With the establishment of the proper adsorption relationship for the freeze-drying simulation,it is possible now to evaluate the effect of initial porosity on the overall drying performance numerically.To this end,several samples with different initial porosity,ε(1-S0),or initial saturation,S0,will be numerically tested.Sample mass was all 1.8 g with the same diameter of 14.8 mm but different heights.Detailed operating conditions are the same as listed in Table 2 except for the ε(1-S0)ranging from 0 to 0.69,corresponding to the S0from 1 to 0.28.The final saturation was taken as 0.03[5].

Fig.5.Variations of drying time with initial porosity.

Fig.5 show s the variations of drying time with different initial porosities.It is expected that the drying time decreases with the increase in the initial porosity,consistent with the experimental results[15].It is interesting to note that there is a threshold of the prebuilt porosity.As it can be seen in Fig.5,when ε(1-S0)is smaller than 0.15,there is insignificant decrease in drying time.Afterwards,the increase in the initial porosity results in an almost linearly decreasing drying time.The existence of the threshold should be a balanced effect between the increased mass transfer and heat transfer.When the pore size is too small,the large mass transfer resistance cannot effectively move the sublimated vapour out of the internal of the bulk.Therefore,the overall performance of the freeze drying would be similar to the solid frozen material.It is believed that this threshold initial porosity is a function of operating temperature,pressure as well as the geometry of the frozen material and the position of heat supply.For low operating temperature,high operating pressure and larger diameter frozen material,the threshold initial porosity is expected to increase.At the other end of the initial porosity,or,the one close to 1,it is postulated that a plateau of drying time should be established at an upper limit of initial porosity because the increase in surface area with the solid material will reach its limit.

For the initial porosity larger than the threshold value and smaller than that of upper limit,drying time is seen decreased with the increase in initial porosity.To keep the same sample mass but with varied initial porosity,apparent surfaceareaor apparent volume of the frozen sample would definitely be different with the sample having a larger initial porosity or a smaller initial saturation giving a larger apparent surface area or a larger apparent volume.This apparent surface area of the currently used porous sample(S0=0.28)is almost as three times as that of the solid one(S0=1.00).This should have improved heat transfer moderately for the sample being dried to absorb more radiation energy.Since the intrinsic porosity,ε,would increase to 0.96 for this type of the material compared with the value of 0.88 for the saturated sample.This implies that the surface area occupancy of the solid substrate,1-ε,would decrease from 0.12 to 0.04,which is one third of the saturated counterpart.This leads to a decreased emissivity,e,in Eq.(9)resulting in a decreased specific surface for radiation heat transfer,and almost offsets the effect of larger apparent area.Consequently,the increased apparent surface area for the initially unsaturated sample has little contribution to the enhancement of radiation heat transfer.Hence,the reduced drying time is not resulted from an increased radiation heat transfer.

It is know n that the conventionally solid frozen material is fully saturated with ice crystals initially,which has a small inherent pore size;w hile the initially porous frozen material is partially filled with crystals,which possesses a relatively large pore size due to pre-built void space.For the former material,drying mechanism is accepted to be the sublimation of free moisture accompanied by interface movement,and followed by the desorption of bound moisture.For the latter one,w hen the initial saturation is large,drying mechanism may be similar to that of conventional freeze-drying.When it is sufficiently small,the pore space occupied by gas has formed continuousn micropore channels during the freezing stage,which leads to a different drying mechanism as show n in previous section although there still exists a sublimation interface.It is because of the pre-built pores pace and larger internal surface area that sublimated and/or desorbed vapour would be easier to escape out of the material being dried,causing an increased drying rate and a decreased mass transfer resistance.That is most likely the primary reason w hy freeze-drying time can be largely saved with the initially porous frozen material.

6.Conclusions

A tw o-dimensional mathematical model of the cylindrical coordinate system was derived to describe the coupled heat and mass transfer be haviours of freeze-drying.Governing equations were solved numerically with the finite deference method.The ADI algorithm can be used to solve the nonlinear equation set with a penta-diagonal coefficient matrix.Numerical simulation con firms that the sublimation interface can be identified naturally;it does not need to be specified as an additional boundary condition in the numerical simulation.

The simple power-law adsorption-desorption equilibrium relation and the Redhead recommended polynomial relation could not represent the experimental drying kinetics curves accurately.The Kelvin's style relation was found to describe the hygroscopic effect satisfactorily during freeze-drying of mannitol aqueous solution,with excellent agreements between the numerical predictions and the experimental data when the parameter,γ,of 5000 was applied.Analyses of temperature and ice saturation pro files show that there is indeed the sublimation-desublimation throughout the entire volume of the initially porous frozen sample.However,the sublimation interface still exists with some residual moisture governed by the adsorption equilibrium in the post-interface region,and additional heat input is required to remove it in the secondary drying stage.Numerical results show that there is a threshold initial porosity for the drying time to decrease significantly with the increase in the initial porosity.The decrease in drying time is attributed to the decreased mass transfer resistance resulted from the increased pore space and internal surface area.

Nomenclature

Subscripts

Superscripts

Acknowledgements

The authors would like to thank Miss Wei Wang sincerely for her help in conducting the calculation work.

[1]M.Brulls,A.Rasmuson,Heat transfer in vial lyophilization,Int.J.Pharm.246(1-2)(2002)1-16.

[2]W.Wang,G.Chen,Numerical investigation on dielectric material assisted microw ave freeze-drying of aqueous mannitol solution,Dry.Technol.21(6)(2003)995-1017.

[3]A.I.Liapis,R.Bruttini,Freeze drying,in:A.S.Mujumdar(Ed.),Handbook of Industrial Drying,second ed.Marcel Dekker,New York 1995,pp.305-343.

[4]Z.H.Wang,M.H.Shi,Numerical study on sublimation-condensation phenomena during microwave freeze drying,Chem.Eng.Sci.53(18)(1998)3189-3197.

[5]W.Wang,G.Chen,Freeze drying with dielectric-material-assisted microwave heating,AICHE J.53(12)(2007)3077-3088.

[6]A.I.Liapis,M.J.Pikal,R.Bruttini,Research and development needs and opportunities in freeze drying,Dry.Technol.14(6)(1996)1265-1300.

[7]R.W.Hartel,Ice crystallization during the manufacture of ice cream,Trends Food Sci.Technol.7(10)(1996)315-321.

[8]A.B.Russel,P.E.Cheney,S.D.Wantling,Influence of freezing conditions on ice crystallization in ice cream,J.Food Eng.39(2)(1999)179-191.

[9]M.J.Millman,A.I.Liapis,J.M.Marchello,An analysis of the lyophilization process using a sorption-sublimation model and various operational policies,AICHE J.31(10)(1985)1594-1604.

[10]J.H.Petropoulos,J.K.Petrou,A.I.Liapis,Network model investigation of gas-transport in bidisperse porous adsorbents,Ind.Eng.Chem.Res.30(6)(1991)1281-1289.

[11]M.J.Pikal,S.Shah,D.Senior,J.E.Lang,Physical chemistry of freeze-drying:measurement of sublimation rates for frozen aqueous solutions by a microbalance technique,J.Pharm.Sci.72(6)(1983)635-650.

[12]S.Nail,L.Gatlin,Freeze drying:principle and practice,in:A.Avis,A.Liebermann,L.Lachmann(Eds.),Pharmaceutical Dosage Forms,Marcel Dekker,New York 1993,pp.163-333.

[13]K.Nakagawa,A.Hottot,S.Vessot,J.Andrieu,Modeling of freezing step during freeze-drying of drugs in vials,AICHE J.53(5)(2007)1362-1372.

[14]W.Wang,Dielectric-material-assisted Microwave Heating in Freeze-drying(Ph.D.Thesis)HKUST,Hong Kong,2005.

[15]K.Yu,W.Wang,Y.Q.Pan,W.Wang,G.Chen,Effect of initially unsaturated porous frozen material on freeze-drying,CIESC J.64(9)(2013)3110-3116(in Chinese).

[16]W.Wang,D.P.Hu,Y.Q.Pan,Y.Q.Zhao,G.Chen,Freeze-drying of aqueous solution frozen with prebuilt pores,AIChE J.61(6)(2015)2048-2057.

[17]W.Wang,G.Chen,Heat and mass transfer in dielectric-material-assisted microwave freeze-drying of skim milk with hygroscopic effect,Chem.Eng.Sci.60(23)(2005)6542-6550.

[18]H.W.Wu,Z.Tao,G.Chen,H.W.Deng,G.Q.Xu,S.T.Ding,Conjugate heat and mass transfer process within porous media with dielectric cores in microwave freeze drying,Chem.Eng.Sci.59(14)(2004)2921-2928.

[19]M.Chen,W.Wang,Y.Q.Pan,K.Yu,G.Chen,1D and 2D numerical verification on freeze drying of initially porous material frozen from aqueous solution,Proceedings of 18th International Drying Symposium(IDS2012),Xiamen,China,2012.

[20]C.J.Geankoplis,Transport Processes and Unit Operations,3rd ed.Engel wood Cliffs,New York,1993.

[21]R.M.Hartel,Phase transitions in ice cream,in:M.A.Rao,R.W.Hartel(Eds.),Phase/State Transition in Foods,Marcel Dekker,New York 1998,pp.327-331.

[22]D.D.Do,Pure Component Adsorption in Microporous Solids,Adsorption Analysis:Equilibria and Kinetics,Series on Chemical Engineering,vol.2,Imperial College Press,London,1998 188-190.

[23]P.Rajniak,R.T.Yang,Unified network model for diffusion of condensable vapors in porous media,AICHE J.42(2)(1996)319-331.

[24]R.B.Bird,W.E.Stewart,E.N.Light foot,Transport Phenomena,second ed.John Willey and Sons,New York,2002 488-508.

[25]E.R.G.Eckert,R.M.Drake,Thermophysical Properties,Analysis of Heat and Mass Transfer,McGraw-Hill,New York,1987 767-783.

[26]P.E.Liley,Physical and chemical data,in:H.R.Perry,D.W.Green,J.O.Maloney(Eds.),Perry's Chemical Engineers' Handbook,seventh ed.Mc Graw-Hill,New York 1997,pp.40-304.

[27]H.Kakiuchi,M.Yamazaki,M.Yabe,S.Chihara,Y.Terunuma,Y.Sakata,T.Usami,A study of erythritol as phase change material,PCMs and Chemical Reactions for Thermal Energy Storage,Second Workshop,So fia,Bulgaria,1998.

[28]I.E.Idelchik,General information and element of aerodynamics and hydraulic of pressure systems,in:M.O.Steinberg,G.R.Malyavskaya,O.G.Martynenko(Eds.),Handbook of Hydraulic Resistance,third ed.CRC Press,Florida 1994,pp.9-12.

[29]D.R.Lide,G.W.A.Milne,Handbook of Data on Organic Compounds,third ed.CRC Press,Florida,1994.

[30]T.Cebeci,J.P.Shao,F.Kafyeke,E.Laurendeau,Computational Fluid Dynamics for Engineers,Horizons Publishing Inc.,Long Beach,2005.

[31]W.Wang,M.Chen,G.Chen,Issues in freeze drying of aqueous solutions,Chin.J.Chem.Eng.20(3)(2012)551-559.