Modelling of industrial processes using computational fluid dynamics

\

Canadian Metallur`ical Quarterly\ Vol[ 26\ No[ 2Ð3\ pp[ 140Ð152\ 0887 Þ 0887 Canadian Institute of Mining and Metallurgy[ Published by Elsevier Science Ltd Printed in Great Britain[ All rights reserved 9997Ð3322:87 ,08[99¦9[99

Pergamon

PII ] S9997Ð3322"86#99926Ð1

MODELLING OF INDUSTRIAL PROCESSES USING COMPUTATIONAL FLUID DYNAMICS MARTHA SALCUDEAN Department of Mechanical Engineering\ University of British Columbia\ 1213 Main Mall\ Vancouver\ B[C[\ Canada\ V5T 0Z3 "Received 08 Au`ust 0885 ^ in revised form 6 July 0886# Abstract*Some examples of computational ~uid dynamics in modelling and optimization of industrial processes are discussed[ Examples include _lm cooling of turbine blades\ gallium arsenide crystal growth\ and black liquor recovery boilers[ Modelling aspects and numerical techniques used are discussed together with some limitations and possible remedies[ Þ 0887 Canadian Institute of Mining and Metallurgy[ Published by Elsevier Science Ltd[ All rights reserved[ Resume*On discute de quelques exemples de calcul par ordinateur de la dynamique des ~uides\ dans la modelisation et l|optimisation de procedes industriels[ Les exemples incluent le refroidissement pelliculaire d|aubes de turbine\ la croissance de cristal d|arseniure de gallium et les fours de recuperation de liqueur noire[ On discute des aspects de la modelisation et des techniques numeriques utilisees ainsi que de certaines limitations et de remedes possibles[ Þ 0887 Canadian Institute of Mining and Metallurgy[ Published by Elsevier Science Ltd[ All rights reserved[

INTRODUCTION In 0835 John von Neumann wrote an enthusiastic assessment of the role computers would play in the future in solving com! plex ~uid ~ow problems[ {{The purpose of the experiments is not to verify a proposed theory but to replace a computation from an unquestioned theory by direct measurement[|| Therefore\ von Neumann hoped that\ with advances in Com! putational Fluid Dynamics "CFD#\ the experiments could be eliminated as the solution of the complex ~ow problems would become possible[ Fifty years after von Neumann made this optimistic assessment of the computer|s role in ~uid dynamics\ the need for experiments has not been eliminated nor is it likely to be in the future[ However\ very signi_cant progress has been made in replacing a large number of experiments\ in guiding experiments and in the expansion of our physical understanding and predictive capabilities[ Researchers are still working very hard at improving computational methods\ and engineers and scientists are working at improving the models of complex processes[ CFD has penetrated into industries where it plays an increasing role not only as a research tool but in the simulation and optimization of industrial processes and products[ The Space and Aircraft industries were an early target for CFD applications\ and other areas such as automotive\ chemical\ and metallurgical industries\ as well as _elds such as hydrology and oceanography followed later ð0Ł[ Because it has enormous potential\ we can look for a rapid increase in the use of CFD methods in industry[ However\ there are still some di.culties in computing three!dimensional large domains\ especially in the presence of a wide range of scales "convergence\ turbulence modelling\ issues of ~ow stability\ and the application of combustion in furnaces and boilers\ etc[#[ 140

There are two distinct but related tasks which have to be undertaken in the modelling of industrial processes[ The _rst task is to establish a model\ namely to analyze the physics and chemistry of the process and determine the governing di}er! ential equations[ In order to establish the governing equations\ it is necessary to have a reasonable understanding of the physi! cal and chemical processes taking place[ For instance\ a judge! ment has to be made as to whether or not radiation would be a factor\ and\ therefore\ if there is a need for a radiation model[ Other issues would then arise*is the radiation model a simple one\ is the gas phase important as far as radiation is concerned and should it be included< Answering these questions requires a good understanding of the physics and chemistry of the process involved and is best handled by engineers[ One of the problems which very often has to be addressed is the modelling of turbulence[ Despite the recent progress in the turbulence simulation\ such as more sophisticated models\ large eddy simulations and direct simulations\ engineers still need to use some simpler models which give a reasonably accurate representation while\ at the same time\ assuring that the com! putations are feasible within reasonable time limits[ Another aspect of the modelling of the physics is establishing the boundary conditions[ The boundary conditions are typically Dirichlet!types when the values of the parameters are known\ von Neumann!types when ~uxes are given\ or mixed types[ Similarly\ there may be information on velocities or pressures out of which boundary conditions can be posed[ Nevertheless\ establishing the computational domain and the boundary con! ditions is by no means a simple task in complex engineering problems[ The computational domains can be extraordinarily large and very complex\ and\ therefore\ intractable through a computation[ Thus\ establishing reasonable boundary con! ditions for a tractable computational _eld becomes a major issue[ Some examples of these problems will be given[

141

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

Once the model has been established\ the second process involves the use of adequate computational techniques with the objective of obtaining reasonable accuracy and performance[ The accuracy of the computations is quite important\ even when the physics and chemistry are not well enough known[ {{Unconverged solutions|| are not acceptable because there is absolutely no guarantee that they are actually close to the real solution[ Local resolution is another important aspect in many engineering problems[ Many large problems contain a very wide range of scales\ and\ therefore\ single grids cannot be used to solve the domain adequately\ as shown in a few examples given later[ It is very important to be able to use local segmentation which allows for a much denser grid in regions in which strong gradients are present[ Another issue which needs to be addressed is the convergence rate of the solution[ It is very important to reduce the computational time for complex prob! lems[ Calculations must converge reasonably\ otherwise\ it is impractical to use CFD as a design and optimization tool[ Using multigrid methods provides very good improvements of the convergence rate for a large number of cases\ an example of which will be presented later[ Turbulence equations con! stitute an on!going challenge in computations[ Turbulence equations typically have very large source terms compared to convective and di}usive terms in many practical situations[ Therefore\ the convergence rate of calculations for turbulence can be poor[ This is particularly true in cases in which body forces are present[ From an engineering point of view and with an awareness of the limitations existing in current turbulence models\ one of the main criteria to be considered when choosing the model is the computational behaviour[ Therefore\ it is very likely that the e}ort in turbulence modelling will address improvements to the model itself in parallel with improving computational e.ciency[ A few examples of process modelling through the use of CFD are given below[ Following the presentation of the problems\ numerical methods are brie~y discussed\ some results are pre! sented\ and a few of the di.culties encountered in the cal! culations are outlined[

FILM COOLING OF TURBINE BLADES Film cooling is a technique used in turbines to protect the blade from the high temperature of the surrounding gas stream[ In this technique\ coolant is injected from a number of discrete\ small injection ori_ces on the blade surface[ The layer of coolant protects the blade surface near the coolant ori_ces and further downstream[ The cooling process involves complex ~ow geo! metries\ high turbulence levels and strong velocity and tem! perature gradients[ A few years ago\ we at UBC undertook both experimental and computational work to investigate the _lm cooling of turbine blades and to optimize cooling geo! metries and other parameters[ The task involved the development of a model which helps to optimize the geometry and injection parameters for the _lm cooling[ The program began a few years ago and has involved both experimental and computational work which has gone through several phases[ An experimental technique has been developed to measure cooling e}ectiveness using a ~ame ion! ization technique and a heat:mass transfer analogy[ Cooling

e}ectiveness {{h|| is de_ned as h  "Taw−T#:"Tc−T# where Taw is the adiabatic wall temperature\ Tc is the coolant tem! perature\ and T is the temperature of the main stream[ Com! putations have been carried out in parallel with experiments for simple geometries\ namely\ for a two!dimensional slot on a ~at plate\ then for a three!dimensional rectangular channel perpendicular to the plate\ and then for channels inclined in di}erent directions[ Finally\ a turbine blade model was mea! sured and calculated[ This step!wise approach is necessary in order to be able too assess the problems associated with the calculations and the experiments ð1Ð00Ł[ Some of the di.culties and problems encountered are pre! sented in the following ] "a# One area of di.culty worth mentioning is the de_nition of the computational domain[ In the real turbine blade\ the computational domain should involve the interior of the blade which contains the internal cooling\ the blade itself through which heat is conducted along with the coolant ori_ces\ and then the blade surface and the main stream[ It is obvious that we are many years away from being able to compute this problem through completely[ Therefore\ this problem has to be separated in such a way that important engineering conclusions can be drawn from the information obtained[ This is not an easy task[ Nearly all investigation done previous to that reported in the present paper started at the blade surface\ assuming a uniform velocity at the exist of the ori_ce[ However\ our calculations and experi! ments have evidenced a strongly non!uniform ~ow at the exit of the coolant ori_ce[ This non!uniformity leads to the conclusion that if the coolant ori_ce is not included in the calculations\ erroneous results are obtained\ mainly closer to the coolant ori_ces which\ in many cases\ are extremely important in leading edge cooling[ Therefore\ the coolant ori_ce was included in the calculation[ However\ further comparison between the experiments and the calculations have shown that the entrance into the ori_ce from the plenum is also important[ Therefore\ part of the plenum should also be included in the coolant ori_ces\ but as the ~ow in the plenum is not well known\ the whole plenum would have to be included[ In consequence\ the size of the domain would become prohibitive[ Therefore\ the coolant ori_ce was included in the computational domain and\ for some cases\ the calculations were expanded into the plenum[ After calculating the simpler two! and three!dimensional planes\ computations were carried out for a UBC exper! imental model turbine blade[ The model has a semi!circular leading edge with four rows of _lm cooling ori_ces posi! tioned symmetrically about the stagnation line as shown in Fig[ 0[ It uses a lateral injection and the cooling ori_ces are inclined spanwise to 29> to the turbine blade surface[ It became very clear from the simple geometries\ that in order to obtain reasonable results\ it is absolutely necessary to resolve the near wall region accurately[ Otherwise\ the sep! aration which takes place at the trailing edge of the ori_ce is missed ð1Ð3Ł[ "b# The numerical method must allow for grid segmentation in such a way that local resolution can be very high[ However\ the consequence of this very _ne mesh is that the standard wall treatment in the kÐo model cannot be used as this

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

142

Fig[ 0[ Illustration of the turbine blade model[ Fig[ 1[ Illustration of staggered grid arrangement[

treatment assumes that the _rst nodal point is located in the fully turbulent region of the ~ow[ This certainly is not the case when a very _ne mesh is used[ Therefore\ the introduction of a low Reynolds number k model near the wall is necessary[ Di}erent turbulence models have been tested\ such as the multiple!time!scale turbulence model[ The results showed some improvements over the standard kÐo model when compared with experimental results[ How! ever\ the most dramatic improvements were due to the increased local resolution in combination with the low Reynolds number k equation ð1Ł[ It is interesting to note that\ at high mass ~ow ratios\ all the improvements failed to assure the necessary accuracy of the modelling[ The results of the calculations at high mass ~ow ratios show signi_cant disagreement with the experiments\ even in the presence of the corrections noted above[ It appears that\ at high mass ~ow ratios\ the ~ow separates and the separation region is not properly resolved by the turbulence models which were tried[ Signi_cant gross unsteadiness\ con_rmed through ~ow visualization\ occurs with enhanced mixing[ This enhanced mixing is not captured by the turbulence modelling[ Another observation noted during the computations was that the jets spread more in the experiments than in the computations[ The reason for this is that the standard tur! bulence kÐo model assumes an isotropic turbulence while the real case turbulence is non!isotropic with the normal component smaller than the spanwise component[ This problem was addressed and the results were somewhat improved by using Bergeles non!isotropic turbulent e}ec! tiveness model ð2Ł[ Another problem area which was encountered was related to the heat transfer calculations in the separation region[ The Reynolds analogy is not valid in that region and\ therefore\ the computations do not adequately rep! resent the real phenomena[ Computational method As shown in Fig[ 1\ a discretized staggered grid arrangement was adopted in which the velocity components are located on the control volume surface and pressures are located at the control volume centre[ The physical tangential velocity com! ponents are used as the primary velocity variables in the momentum equation[ The numerical scheme is formulated directly using the co!ordinate invariant governing equations\ i[e[\ the equations are expressed with operators\ gradient D and

divergences 9\ without referring to a co!ordinate system[ For any given grid\ the physical control volumes are constituted by using the discrete grid nodes as their vertices[ The co!ordinate invariant governing equations are then integrated directly in the physical control volumes[ The discrete governing equations are obtained by expressing the resulting surface ~uxes using the discrete velocity values and the grid physical geometric quantities\ which include the volumes\ the surface areas and the surface normal directions of the physical control volumes[ These physical geometric quantities are calculated in such a way that they satisfy the geometric conservation laws[ The conservation equations used are the following ] ru = 9u−9 = "me}9u#  −9p\

"0#

9u  9\

"1#

mt 9k  G−ro sk

1

"2#

mt o o1 9o  C0 G−rC1 \ so k k

"3#

0

9 = ruk−

0

9 = ruo−

0

9 = ruT−

1

0

1 1

mt ml 9T  9\ ¦ Pr Prt

"4#

where k is the turbulence kinetic energy\ o is the turbulence energy dissipation rate\ r is the ~uid density\ G is the turbulence energy generation rate\ T is temperature\ Pr and Prt are the Prandtl number and turbulent Prandtl number\ respectively\ and C0\ C1\ sk and so are empirical constants as recommended by Launder and Spalding "0863#[ A domain segmentation technique is used[ The method divides the domain of interest into di}erent ~ow regions[ The ~ow domain over a turbine blade is decomposed into the main ~ow region and a number of injection tube regions[ In the present system\ attention is focussed on developing a scheme to allow the use of non!smooth grids with rapidly changing slopes across the interface[ This was achieved by manipulating the discretized momentum equations at the interface to assure the momentum conservation across it[ The method is particularly useful in dealing with inclined injection ori_ces ð4Ð6Ł[ A multigrid method has been developed and used to accel! erate the convergence rates[ There are additional di.culties associated with the multigrid calculation of turbulent ~ow in curvilinear grids[ A mass!conserving restriction is used for the

143

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

Fig[ 2[ Illustration of the turbine blade model used in the experimental study[

transfer of the velocity components\ while the full weighting is applied for the restriction of residuals[ A bi!linear interpolation is adopted as a prolongation operator[ Based on the method outlined\ a three!dimensional CFD code\ CMGFD\ was developed using the multigrid method in multi!zone curvilinear grids[ This code has the capacity to solve incompressible tur! bulent ~ows in arbitrary three!dimensional geometries once the curvilinear grids are generated[ An elliptic grid generation method proposed by Thompson is used to generate the grid[ In order to obtain an e.cient grid generator\ a multigrid method was developed to solve these equations[ The Full Approximation Scheme "FAS#\ in con! junction with the Full Multigrid Cycling Scheme "FMG# is adopted[ The basic idea of the FMG Scheme is to provide a good initial solution for the _ne grids by prolongating from the solution on coarse grids[ The multigrid iteration procedure is carried out using the V!cycle[ The turbine blade model used in the experimental study is illustrated in Fig[ 2[ The computational domain includes the main ~ow region\ with two injection ori_ces in the _rst row and one injection ori_ce in the second\ as illustrated in Fig[ 3[ The

Fig[ 4[ Coolant main stream mass ~ow ration M  9[41 ^ solutions at the streamwise plane "xz!plane# through the centreline of the 04> coolant ori_ce ] "a# local velocity _led around the curved blade surface ^ "b# local velocity _eld near the ori_ce exit[

main ~ow region includes 49d upstream from the leading edge\ 39d downstream in the positive X direction from the end of the semi!cylinder and 29d in the Z axis direction from the leading edge[ The _ne grid used for the study contains a 090×06×26 nodes and 8×8×08 nodes for each injection duct[ The periodic "cyclic# condition is used on the boundaries in the spanwise direction[ It was observed that heavy under!relaxation was necessary in order to obtain a converged solution ð6Ł[ Some of the results are illustrated in Figs 4"a# and 4"b#\ 5 and 6[

Fig[ 3[ Illustration of the computational domain for the staggered double row injection\ including a main ~ow region\ two half injection!hole ducts in the _rst row and one injection hole duct in the second row[

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

144

Fig[ 5[ Contours of predicted _lm cooling e}ectiveness for coolant main stream mass ~ow ratio M  9[41 at the blade surface[

experiments were carried out to assess the e}ect of density di}erences ð6Ł[ Carbon dioxide was used as a coolant in both experiments and computations[ The e}ect of the density variation on the cooling e}ectiveness has been assessed[ Based on the experiments and computations carried out\ an optimum geometry with excellent cooling e}ectiveness characteristics has been identi_ed and is being implemented by Pratt + Whitney Canada[

MODELLING OF MAGNETIC CZOCHRALSKI GROWTH OF SEMI!CONDUCTOR CRYSTALS Fig[ 6[ M  9[41 ^ comparison of the spanwise averaged _lm cooling e}ectiveness with the measured result[

Conclusions based on these _gures and on others which can! not be included due to space constraints\ are outlined in the following ] "0# The ~ow _elds show signi_cant di}erences between the _rst and second rows because of the pressure variation on the blade surface[ "1# Flow at the exit of the ori_ces is very non!uniform\ par! ticularly in the front row[ "2# The interaction between the main ~ow and coolant shows a di}erent vortex formation in the spanwise injection from that observed in streamwise injection[ "3# The vortices are weaker when formed further downstream and are further from the blade surface\ and therefore do not produce the lifting of the coolant from the surface which occurs for streamwise injection[ "4# Examination of the _lm cooling e}ectiveness shows a good coverage at coolant main stream mass ~ow ratios of M  nc:n  9[41 "M being the mass ~ow ratio averaged over all coolant holes for double row cooling# where the coolant from the _rst row ~ows between the ori_ces of the second row[ For M  9[86\ the coolant from the _rst row of ori_ces blows towards the ori_ces of the second row and the coverage is therefore poor[ A comparison of _lm cooling e}ectiveness with experimental observations shows a very good agreement at M  9[41 and somewhat deteriorating agreement for M  9[86[ Further computations and

Gallium Arsenide "GaAs# crystals needed for micro!elec! tronic applications are mainly produced by the liquid encap! sulated Czochralski technique "LEC#[ This technique "Fig[ 7# involves pulling a single crystal slowly from molten Gallium Arsenide through a layer of liquid encapsulant and into a cooler stainless steel vessel containing argon at a high pressure "04Ð19 atm#[ The melt and encapsulant\ contained in a pyrolytic boron nitride crucible\ are heated from the side by resistance heaters[ The encapsulant\ usually boric oxide\ prevents the loss of vola!

Fig[ 7[ Schematic of LEC growth of GaAs[

145

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

Fig[ 8[ Flow and temperature _eld in the puller for a crystal rotation rate of 19 rpm and a crucible counter! rotation rate of 5 rpm for a melt height of 9[94 m when no magnetic _eld is applied[

tile arsenic from the melt and from the hotter surface of the crystal near the crystalÐmelt interface[ The crystal and crucible are rotated to minimize the thermal asymmetries in the system\ and the crucible is raised to maintain the crystalÐmelt interface at the speci_ed region of the heater[ LEC!grown GaAs crystals often have a large number of dislocations and striations of impurities[ The electrical charac! teristics of GaAs devices are critically a}ected by the level and distribution of these defects in the crystals[ Temperature ~uctuations in the melt arising from ~ow oscillations often cause the striations of impurities in the crystal[ Because molten GaAs has a large electrical conductivity\ the thermal convection in the melt can be suppressed by the application of a magnetic _eld[ Flow suppression also helps to control the temperature gradients near the crystalÐmelt interface and to maintain a nearly ~at interface[ However\ ~ow near the crystalÐmelt inter! face is necessary for the impurities segregated at the interface to be transported to the bulk of the melt[ Therefore\ knowledge and understanding of thermal processes in the LEC puller with and without the application of a magnetic _eld are essential for growing better quality crystals[ Mathematical model Flow and heat transfer processes in the puller during LEC growth of GaAs are extremely complex[ In the following\ the ~ow and heat transfer in the melt and encapsulant\ and heat transfer in the crystal are considered[ A simple model for ther! mal radiation in the puller is used[ The melt\ encapsulant and crystal are assumed to be gray bodies\ and each of them is assumed to exchange radiation only with an ambient environ! ment at a constant and steady temperature\ T[ The ~ow and heat transfer in the ambient gas are neglected[ A steady\ uniform axial magnetic _eld\ B9\ is assumed to be applied for ~ow suppression in the melt[ The e}ects of an induced magnetic _eld on the applied _eld are neglected because the magnetic Reynolds number is much smaller than 0 "less than 09−2# for the conditions encountered during the LEC

growth of GaAs[ The crystal is considered to be electrically! conductive[ The crystal diameter is kept constant\ the crucible side wall is assumed to be at a uniform and constant temperature\ Tc\ and the pull rate is determined from an energy balance at the crystalÐmelt interface[ The crucible bottom wall is assumed to be insulated[ The top of the crystal is assumed to be at the ambient temperature\ T[ The defects in the crystal are lower when the growth con! ditions are adjusted to maintain a nearly ~at crystalÐmelt inter! face[ Therefore\ in this study\ the crystalÐmelt and the melt! encapsulant interfaces are assumed to be ~at and at the same height[ The thermocapillary e}ects due to the gradients in the melt!encapsulant interfacial tension are assumed to be negli! gible[ The system is assumed to be axisymmetric[ Further details of the mathematical model\ governing equations and boundary conditions\ and the numerical method can be found in ð01\ 02Ł[ Results and discussion In LEC growth of GaAs\ the crystal and crucible are normally counter!rotated to reduce the thermal asymmetries in the system and to form a boundary layer of nearly uniform thickness in the melt beneath the crystal[ In MLEC growth of GaAs\ the e}ects of crystal and crucible rotations are not clear[ We obtained the ~ow and heat transfer in the melt and encapsulant\ and heat transfer in the crystal with and without the application of an axial magnetic _eld for a crystal diameter of 9[093 m\ a crucible diameter of 9[073 m and melt heights of 9[94 and 9[915 m[ Various combinations of crystal and crucible rotational rates from 9 to 19 rpm were considered[ The e}ects of both co! rotation and counter!rotation of the crystal and crucible were examined[ For a typical combination of vs  19 rpm and vc  −5 rpm\ and for a melt height of 9[94 m\ the ~ow and temperature _elds with and without the application of an axial magnetic _eld are shown in Figs 09 and 00[ When the magnetic _eld is not applied\ the meridional ~ow in the melt is bicellular[ The smaller cell near

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

146

Fig[ 09[ Flow and temperature _eld in the puller for a crystal rotation rate of 19 rpm and a crucible counter! rotation rate of 5 rpm for a melt height of 9[94 m with an axial magnetic _eld of 9[4T[

Fig[ 00[ Flow and temperature _eld in the puller for a crystal rotation rate of 19 rpm and a crucible counter! rotation rate of 5 rpm for a melt height of 9[915 m when no magnetic _eld is applied[

the crystalÐmelt interface occurs due to the crystal rotation[ The larger cell occupying nearly the entire melt occurs due to the crucible rotation[ Though vs is much larger than vc\ the e}ects of crystal rotation are con_ned to only a small region beneath the crystal[ The nearly ~at isotherms beneath the crystal indicate that the crystalÐmelt interface is likely to be ~at for high crystal rotational rates[ However\ the small radial temperature gradi! ents near the crystal\ melt and encapsulant tri!junction indicate that the diameter control is likely to be di.cult[ Due to ~ow suppression "Fig[ 09#\ the heat transfer in the melt is dominated more by conduction when a magnetic _eld is applied[ The crucible side wall temperature needs to be increased with the strength of the magnetic _eld to avoid freez! ing the top surface of the melt[ Though the meridional ~ow in the melt is substantially suppressed\ the e}ects of convection are still important[ As a result\ the value of Tc for satisfactory crystal growth for B9  9[4T is lower than that predicted by the conduction analysis[ The shapes of isotherms in the melt near the crystalÐmelt interface indicate that the interface is likely to be ~atter and the diameter control is likely to be easier with the application of a magnetic _eld[ The ~ow in the encapsulant is unicellular[ The encapsulant is

an electrical insulator\ and\ therefore\ the meridional ~ow in it is not suppressed by the magnetic _eld[ The meridional ~ow circulation strength in the encapsulant is higher in Fig[ 09 than in Fig[ 8 due to a larger value of Tc for B9  9[4T[ As the melt height decreases due to the growth of the crystal\ the crucible is normally raised so that the crystalÐmelt interface is in a speci_ed region of the puller[ As a result\ for lower melt heights\ the crucible bottom is more exposed to the radiant heat from the heater[ Therefore\ for the calculations with lower melt heights\ the bottom of the crucible is assumed to be at the same temperature as that of the crucible side[ Figures 00 and 01 show the ~ow and temperature _elds with and without the application of a magnetic _eld for a melt height of 9[914 m[ When no magnetic _eld is applied\ the calculations show that the ~ow and temperature _elds in the melt are strongly dependent on crystal and crucible rotational rates[ For the combination of the rotational rates shown in Fig[ 00\ the e}ects of crystal rotation on the ~ow in the melt are more pronounced as seen from the size of the cell beneath the crystal[ The iso! therms indicate that the diameter control will be di.cult when no magnetic _eld is applied[ The ~ow in the melt is suppressed by the magnetic _eld and is less sensitive to the rotational rates

147

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

Fig[ 01[ Flow and temperature _eld in the puller for a crystal rotation rate of 19 rpm and a crucible counter! rotation rate of 5 rpm for a melt height of 9[915 m with an axial magnetic _eld of 9[4T[

when an axial magnetic _eld is applied[ The ~ow\ however\ becomes multi!cellular when the magnetic _eld is applied[ The shapes of isotherms in the melt re~ect the nature of multi! cellular ~ow\ and indicate that the crystalÐmelt interface is {{wavy|| and that the diameter control may be more di.cult with the application of a magnetic _eld[ This shows that an optimization of the rotational rates and the magnetic _eld strength is needed to obtain ~at crystalÐmelt interface with better diameter control[ Of several problems encountered during the modelling of the crystal growth\ some were related to the very large domain of calculation which basically includes the liquid Gallium Arsen! ide\ the solid crystal\ the encapsulant and the gas in the presence of several heat transfer mechanisms\ namely\ forced and natural convection in the liquid metal\ electro!magnetic convection\ occasionally surface tension induced convection\ and natural convection in the gas[ Another major di.culty is related to calculating the interface\ which requires domain segmentation as well as curvilinear co!ordinates[ An interesting observation was made that\ under certain growth conditions\ the ~ow became unstable and bifurcation occurred[ The ~ow exhibited solutions with a varying number of cells ^ for instance\ di}erent solutions were noticed for the same conditions exhibiting a di}erent number of cells[ A transient solution shows a con! tinuous variation ^ this is very important as it indicates unstable conditions for the crystal growth and in~uence the crystal quality signi_cantly[

FLOW IN RECOVERY BOILERS Kraft recovery boilers are large "typically 09×09×29 m# and expensive installations used in pulp mills to generate heat from black liquor and to convert inorganic chemicals in the liquor to a form suitable for re!use in the pulping process[ Black liquor is introduced through liquor guns while the air is introduced through ports[ Typical recovery boiler designs have two or three levels of air ports\ the lower level being called primary air ports\ the mid!level\ secondary air ports and the upper level\ tertiary air ports[ The capacity of recovery boilers is e}ected mainly by

carry!over and fume!related plugging[ Fume liquor droplets and char particulates become entrained in the high velocity gas ~ows which causes the plugging of screens\ tubes\ superheaters and the generating bank[ Therefore\ boilers sometimes have to be operated at reduced loads[ Optimizing such units can increase their capacity\ thermal output\ and running time between water washes[ Recovery boiler performance is critically dependent on the air and liquor delivery systems[ Many di}er! ent options are promoted for these systems by boiler suppliers[ There is little data to support these options as\ until recently\ their design was mainly based on intuition and experience[ Knowledge of the ~ow and combustion in recovery boilers is very di.cult to obtain[ Detailed measurements in an operating boiler are prohibitive because of the hostile environment[ Even cold ~ow measurements are di.cult to obtain due to the large size of the unit and are expensive due to the cost associated with the loss of production during the measurement work[ Mathematical modelling of the ~ow and combustion in a boiler is an attractive alternative as the required details can be obtained and the e}ects of the operating parameters can be investigated without in!boiler measurements[ The recovery boiler model development program has been underway since 0889 at the University of British Columbia[ At present\ the hot ~ow model has been fully developed and is at the commercial application stage[ The modelling has involved several tasks[ One was the development of mathematical models to simulate the hot ~ow combustion and chemistry as well as measurements\ both in actual recovery boilers and in physical models of boilers[ Both mathematical and physical modelling were carried out in an integrated manner\ and full!scale measurements were carried out by Paprican in close collabor! ation with UBC ð22Ł[ The physical modelling involves 0] 19 scale plexiglass models in which water is used as a working ~uid to provide visualization of the ~ow[ Flow visualization studies have been carried out using a light sheet technique to determine the overall ~ow patterns\ asymmetries and instabilities[ Detailed ~ow velocity measurements and turbulence characteristics have been obtained using a Laser Doppler Velocimeter "LDV#[ Com! parisons have been carried out between the water ~ow measure! ments and computations ð09\ 17\ 23\ 25Ł[

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

Mathematical modellin` In mathematical modelling\ the ~ow _elds were obtained by solving the governing equations representing the conservation of mass\ momentum and energy of the gas[ The ~ow in recovery boilers is turbulent and\ with variable density\ conservation equations also have to be solved to simulate the e}ects of turbulence and to determine the temperature\ and\ therefore\ a radiation model is needed[ The gas species concentration and the liquor product particulate distribution also has to be deter! mined ð18\ 29Ł[ The output of the model includes descriptions of ] "i# the gas ~ow _elds within the boiler cavity\ including the velocity _eld and concentration of O1\ H1O and CO1\ and combustibles "H1CO and CH3# ^ "ii# particle burning\ including distribution of states "e[g[\ drying\ pyrolysis\ char and smelt# in each region of the boiler ^ "iii# the amount of carry!over of burning particles and en! trained smelt "inorganic material# carried out of the boiler cavity with the combustion gases ^ and "iv# temperature distribution[ The mathematical modelling of ~ows in recovery boilers is a complex task[ The ~ow is fully three!dimensional as geometric symmetry\ even on one plane\ is not present in all designs[ Even when the boiler is geometrically symmetric\ ~ow instabilities and jet interactions can cause major asymmetries in the ~ow patterns[ The ~ow is often grossly unsteady and such unsteadi! ness has been observed in physical models of boilers as well[ The turbulence is treated using the kÐo model[ A log normal distribution is assumed for the liquor droplet diameters[ The black liquor spray is simulated using a method proposed by Frederick and Grace ð03\ 04Ł[ Each droplet is tracked from the time of injection until it lands on the bed or is carried over[ The exchange of mass and energy between the droplet phase and the gas phase is taken into account[ The gas combustion rate is computed using the Magnussen model ð05Ł[ The gas radiation heat transfer is calculated using a model developed based on the discrete!ordinates method ð06Ł[ The system of equations solved includes the conservation equations for the mass\ momentum\ energy and species "O1\ CO\ CO1\ H1O\ H1\ CH3\ N1# concentration[ In addition\ the conservation equations for the turbulence quantities "k and o# and the equations for the radiative intensities are solved[ The discrete ordinates method is used for the radiation heat transfer computations[ The radiation model takes into account the emit! ting\ absorbing and scattering medium such as the particle! laden gas present in a recovery boiler[ The exponential wide! band model is used for the prediction of radiative properties of combustion gases\ and Mie scattering is used for the prediction of radiative properties of particles[ In a recovery boiler there are millions of black liquor droplets in suspension at any given moment[ A statistical method is used to simulate these droplets in the computer model[ Typically\ 09\999 to 49\999 computational droplets are used[ The results are independent of the number of computational droplets[ Each computational droplet is assumed to represent a number of real droplets that have the same location\ velocity\ chemical composition and physical properties[ A log normal distribution

148

is assumed for the droplet diameters at injection[ The mean droplet diameter for this distribution is taken from measure! ments carried out at IPST ð07Ł[ The momentum equation for each droplet is solved\ taking into account the drag\ weight and inertial forces[ Each droplet undergoes drying\ pyrolysis\ char burning and smelt reduction phases[ The droplet releases water to the gas phase during the drying stage[ During pyrolysis\ the droplet releases combustibles "CO\ CH3\ etc[# to the gas phase[ Both drying and pyrolysis are assumed to be heat!transfer con! trolled[ The char carbon gasi_cation "burning# is assumed to occur via reactions with O1 and H1O[ The char burning rate is assumed to be controlled by the rate of mass transfer of O1 and H1O to the char particle[ The volatile combustion model includes burning of CH3\ CO and H1[ The combustion is assumed to be controlled by the mixing rate\ as the typical temperature in a recovery boiler is higher than 0999>C "the Arrhenius expression does not apply#[ The Magnussen ð05Ł com! bustion model is used to compute the reaction rates as a func! tion of the turbulent kinetic energy and the turbulent energy dissipation rates[ The black liquor droplet computations are carried out in a moving frame of reference and the gas phase computations are carried out in a _xed Eulerian frame[ The gas phase and droplets exchange mass\ momentum and energy via source sink terms in the appropriate equation[ For example\ if a droplet evaporated water while crossing a cell\ the water vapour mass is added to the cell mass balance through the source term[ Computational procedures The governing di}erential equations described above are inte! grated over the control volumes and\ after application of Gauss| divergence theorem\ integral equations for the convective and di}usive ~uxes across the control volume faces are obtained ð14Ł[ Two major issues critical to the simulation of transport phenomena in recovery boilers were addressed[ The _rst is a need for high local resolution\ especially at the air port levels\ and the second is the need to address slow convergence rates of the solution algorithm for large domains typical of recovery boilers[ The slow convergence of the solution algorithms is due generally to the persistence of low!frequency errors that are not e}ectively removed at a grid which is _ne relative to the wavelengths of the errors[ In order to remove these low! frequency components of the errors\ a multigrid solution algo! rithm was used[ In this technique\ relaxation techniques are applied on a hierarchy of grids\ so that error components cor! responding to a wide range of frequencies are e}ectively removed[ Convergence is obtained rapidly using the multigrid algorithm described in Ref[ ð14Ł[ The local resolution was addressed by using a block!struc! tured segmented grid[ The grid in each segment is completely independent of the grid at the other segments\ and is based on the requirements for the local ~ow resolution and the air port geometry[ Since the grid in each segment is independent\ the control volume faces at the segment interfaces do not necessarily match[ The segmentation feature allows for a representation of the details of the ~ow in critical zones of the recovery boiler\ for example\ where jet interaction takes place at the air port level[ It also allows for a more accurate representation of the air ports\ as the grid of each segment is designed speci_cally to

159

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

take into account the air port sizes and locations in that segment[ The grid distribution is optimized using this technique\ as over 79) of the grid nodes can be concentrated in regions of interest[ This capability is critical as a realistic representation of the jet interaction is needed to provide information for the computation of the droplet trajectories and the chemistry[ It should be noted that the complexity of the computer code increases signi_cantly to provide such ~exibility in de_ning the grid[ The solution method has been developed from a technique originally proposed by Vanka ð19Ł[ This technique has been extensively modi_ed to make rapid computation of recovery boiler ~ows possible[ The modi_cations include additional terms in the momentum equations such that overall mass con! servation is satis_ed\ even in the initial iterations[ In addition\ several terms in the momentum and turbulence equations are damped to prevent divergence ð27Ł[ The modelling and computations of black liquor recovery boilers have raised a signi_cant number of issues[ Some of the modelling issues include the very large number of primary air ports which makes representation of the individual jets very di.cult\ the need to impose a velocity distribution on the air port rather than computing through the windbox\ and the insta! bilities generated by the interaction of jets which have high momentum ratios to the cross ~ow[ Obviously\ these are only some of the di.culties which have been addressed[ In numerical simulations\ closely!spaced primary air ports

are often represented by slots in order to reduce the com! putational requirement[ The slots allow for a coarser grid to be used to represent the jet outlets which are numerous[ This attempt at reducing the computational requirement is very important\ especially in comprehensive furnace simulations where many equations representing various physical and chemi! cal processes need to be solved in each grid cell[ To limit the computational requirement to a manageable level\ the rep! resentation of rows of primary jets by equivalent slot jets is essential[ Simulations have been carried out for simpli_ed geo! metries\ including independent jets[ The results have shown that changing the boundary data for k and o for slot jets can a}ect the resulting predictions[ By properly choosing boundary data for these turbulence quantities\ slot representation can generate ~ow _elds that adequately re~ect the characteristics of the ~ow _eld due to discrete primary jets[ A suggestion for boundary data for use in the slot representation is to choose values for turbulence intensity\ I\ and length scale\ Ld\ to be\ respectively\ about 09 and 20[5 "092:1# times larger than those normally used for a single turbulent jet which constitutes the row of jets ð17Ł[ The complexity of the calculation creates the need to impose a velocity distribution on the air port rather than calculating through the windboxes[ The investigation for the simple wind! box design has shown that\ although the ~ow inside the windbox unit can be very complicated\ nevertheless\ the velocity pro_le at windbox outlet is close to uniform\ and with values that

Fig[ 02[ Generic recovery boiler schematic and air systems[

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

deviate little from the bulk average value computed based on mass conservation[ This observation lends support to the use of simple uniform velocity pro_les in the setting up of boundary conditions for air jets in numerical recovery furnace simulations ð17\ 22Ł[ Instabilities were observed which were generated by the inter! action of jets with high momentum ratios to the cross ~ow[ It was observed that when opposing jets interact\ unsteady turbulent ~ow behaviour occurs[ There is a physical instability when opposing jets collide\ either in the presence of a cross! stream or a closed!~oor cavity[ This instability produces a bifur! cation to several steady state ~ows or an oscillatory ~ow that can be captured by a suitable numerical method[ The cal! culations show a perfectly periodic result for the mean ~ow[ Experiments to con_rm the computations showed a very similar pattern\ with some variability in the oscillations[ This _nding is very important as this ~ow behaviour can lead to unstable and unpredictable ~ow patterns in the boiler\ and\ therefore\ make it impossible to control the process ð29\ 21Ł[ There is a very complex coupling between chemistry\ ~ow and radiation[ Tests with the radiation model have shown that a relatively large number of directions in the radiation model is needed[ The model is much more tolerant of a coarse grid than of a low number of directions[ Because of computational memory and speed\ it was necessary to carry out the radiation calculation on a grid coarser than that used for the hyd! rodynamic calculations[ However\ the number of directions was larger which led to a far better closure of the radiation equation[ Cold ~ow is calculated _rst and\ based on the results\ the droplet trajectory and heat release is then calculated[ Based on this

150

data\ a new ~ow _eld is calculated and the process is iterated until convergence is reached[ The heat sources are very large and sometimes cause instabilities[ The treatment of combustion near the bed is very complex and good resolution is necessary[ There is a lack of understanding of the radiation properties of the particles and of deposition in the heat exchanger region[ The heat exchanger region is very complex and\ therefore\ a simpli_ed treatment based on the assumption of porous medium is necessary[ The model developed can be used as a tool to optimize both air systems and liquor delivery systems[ As an example\ a gen! eric recovery boiler was investigated[ The boiler has a 09 m square base and is 39 m high\ with the bullnose located 15[4 m from the ~oor "Fig[ 02#[ The boiler features three levels of air ] four!wall primary at 0[1 m above the ~oor ^ four!wall secondary at 2 m above the ~oor ^ and two!wall tertiary at 09 m above the ~oor[ The boiler operates with an air split of 37) primary\ 21) secondary and 19) tertiary[ The three air systems inves! tigated were ð15Ł ] "0# a base case ^ "1# the two!wall secondary ^ and "2# interlaced tertiary[ All of these systems are represented in Fig[ 02[ The jet locations are indicated by arrows in all three cases ^ the size and shape of the air ports and the total air injected at each level are the same[ The base case has _ve secondary ports on each wall and seven tertiary ports on the front and rear walls[ In the two! wall secondary case\ the secondary level con_guration has been

Fig[ 03[ Isotherms in the frontÐrear vertical central plane ] "a# base case ^ "b# two wall secondary ^ "c# interlaced tertiary[

151

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES

modi_ed[ The side wall jet velocities have been increased and the front and rear wall jet velocities have been decreased\ respec! tively\ by 13) so that the total air input remains unchanged[ The primary and tertiary levels are the same as in the base case[ In the interlaced tertiary case\ the tertiary level con_guration has been changed[ The velocities of the opposing jet pairs have increased or decreased by 27) to create interlaced tertiary patterns[ The primary and secondary level jets are the same as in the two other cases[ Figure 03 shows the vertical velocity contours in the frontÐrearÐcentre plane for the three di}erent air systems[ Each case shows central upward ~ow and down! ward ~ow close to the walls[ The central upward ~ow core is narrower in the base case and is of higher velocity in the two! wall secondary case[ For the interlaced tertiary case\ the upward ~ow region is pushed farther away from the bullnose\ resulting in a large recirculation region under the bullnose[ The gas temperature is higher in the centre\ lower near the walls and higher just above the bed\ as shown in Fig[ 04[ In order to compare the performance of the air systems\ the carry!over has been computed[ Most of the particles consist of smelt ^ however\ some of the larger particles are still in the char! burning stage of combustion[ The base case has the highest total carry!over of about 01) of the total black liquor mass ~ow[ The carry!over distribution of the two!wall secondary case is somewhat lower and that of the interlaced case is the best[ It can be observed that modifying the secondary and tertiary air input or the spray characteristics can produce changes in the mechanical carry!over[ Evaluations have also been carried out for liquor spray characteristics which are not presented here due to space constraints[

The model has been applied extensively to a large number of industrial boilers for evaluation of di}erent designs and parametric functions[ It provides important information for a more objective evaluation of the boiler characteristics[ Exten! sive ~ow visualization carried out with the computations has proven to be very helpful in understanding the ~ow and has produced a very positive response from the boiler operators[

CONCLUSION Computational Fluid Dynamics has been established as a major tool in process simulation and optimization[ There are major challenges in computing real engineering processes with respect to setting up reasonable physical models as well as in carrying out the computations[ In order for an undertaking to be successful\ the modeller has to have a good understanding of the physical processes[ The computations are very demanding and require skilled personnel[ The non!linearity of the equa! tions\ the complexity\ the three!dimensionality\ the turbulence and other associated transport phenomena of heat and species\ make for a very complex problem which requires a sophisticated analyst[ However\ the costs of highly!quali_ed personnel\ com! puting\ etc[ are certainly minor compared to the very large cost entailed in a trial and error approach to process simulation and optimization[ Therefore\ CFD should be used extensively by process engineers[ Acknowled`ements*This paper is the result of work carried out by a

Fig[ 04[ Upward velocity contours in the front!rear vertical central plane ] "a# base case ^ "b# two wall secondary ^ "c# interlaced tertiary[

M[ SALCUDEAN ] MODELLING OF INDUSTRIAL PROCESSES large group of researchers\ and the author would like to acknowledge the following principal contributors ] Professor Ian Gartshore\ Dr Zia Abdullah\ Dr Paul Nowak\ Dr Peri Sabhapathy\ Paul Matys and Dr Pingfan He[ Funding from the Department of Energy "US#\ Energy Mines and Resources Canada\ Industry Canada\ Pratt + Whitney Canada and NSERC is also gratefully acknowledged[

10[ 11[ 12[

REFERENCES 0[ Salcudean\ M[\ Computational ~uid ~ow and heat transfer*an engineering tool[ Transactions CSME\ 0880\ 04"1#[ 1[ Zhou\ J[!M[\ Salcudean\ M[ and Gartshore\ I[ S[\ Application of the multiple!time!scale turbulence model to the prediction of 1!D _lm cooling e}ectiveness[ Paper presented at the 09th International Heat Transfer Conference\ Brighton\ U[K[\ 03Ð07 August 0883[ 2[ Zhou\ J[!M[\ Salcudean\ M[ and Gartshore\ I[ S[\ Prediction of _lm cooling by discrete!hole injection[ Paper presented at the Inter! national Gas Turbine and Aeroengine Congress and Exposition\ Cincinnati\ Ohio\ 13Ð16 May 0882[ 3[ Zhou\ J[!M[\ Salcudean\ M[ and Gartshore\ I[ S[\ A numerical computation of _lm cooling e}ectiveness[ Near Wall Turbulent Flows[ Elsevier Science Publishers\ 0882[ 4[ He\ P[ and Salcudean\ M[\ A numerical method for 2!D viscous incompressible ~ows using non!orthogonal grids[ International Journal for Numerical Methods in Fluids\ 0883\ 07\ 338Ð358[ 5[ He\ P[\ Salcudean\ M[ and Gartshore\ I[ S[\ Computations of _lm cooling for the leading edge region of a turbine blade model[ Paper presented at the ASME Turbo Expo |84\ Houston\ Texas\ 4Ð7 June 0884[ 6[ He\ P[\ Licu\ D[\ Salcudean\ M[ and Gartshore\ I[ S[\ Leading edge _lm cooling ] computations and experiments including density e}ects[ Paper submitted for presentation at the 0885 ASME Inter! national Gas Turbine and Aeroengine Congress\ Birmingham\ U[K[\ 09Ð02 June\ 0885[ 7[ Salcudean\ M[\ Gartshore\ I[ S[\ Zhang\ K[ and McLean\ I[\ An experimental study of _lm cooling e}ectiveness near the leading edge of a turbine blade[ Journal of Turbomachinery\ 0883\ 005\ 60Ð 68[ 8[ Findlay\ M[\ He\ P[\ Salcudean\ M[ and Gartshore\ I[ S[\ A row of streamwise!inclined jets in cross~ow ] measurements and calcu! lations[ Paper submitted for presentation at the 0885 ASME Inter! national Gas Turbine and Aeroengine Congress\ Birmingham\ U[K[\ 09Ð02 June\ 0885[ 09[ Ajersch\ P[\ Ketler\ S[\ Zhou\ J[!M[\ Gartshore\ I[ S[ and Salcudean\ M[\ Multiple jets in a cross~ow ] detailed measurements and numeri! cal simulations[ Paper presented at the ASME Turbo Expo |84\ Houston\ Texas\ 4Ð7 June 0884[ 00[ Gartshore\ I[\ Licu\ D[ and Salcudean\ M[ E[\ Film cooling for gas turbine blades[ Keynote Paper for the 01th Australasian Fluid Mechanics Conference\ Sydney\ Australia\ 09Ð04 December 0884[ 01[ Sabhapathy\ P[ and Salcudean\ M[ E[\ Journal of Crystal Growth\ 0878\ 86\ 014[ 02[ Sabhapathy\ P[ and Salcudean\ M[ E[\ Journal of Crystal Growth\ 0889\ 093\ 260[ 03[ Frederick\ W[ J[\ Noopila\ T[ and Hupa\ M[\ Combustion behavior of black liquor at high solids _ring[ TAPPI J[\ 0880\ 63"01#[ 04[ Grace\ T[ M[\ Personal communications[ 05[ Magnussen\ B[ F[ and Hjertager\ B[ H[\ On mathematical modeling of turbulent combustion\ 05th Symp[ "Int[# on Combustion[ The Combustion Institute\ Pittsburgh\ 0865\ pp[ 608Ð618[ 06[ Fiveland\ W[ A[\ Three!dimensional radiative heat transfer solu! tions by the discrete!ordinates method[ Journal of Thermophysics and Heat Transfer\ 0877\ 1\ J298ÐJ205[ 07[ Empie\ H[ J[\ Lien\ S[ J[\ Yang\ W[ and Samuels\ D[ B[\ E}ect of black liquor type on droplet formation from commercial spray nozzles[ Journal of Pulp and Paper Science\ 0884\ 10"1#[ 08[ Adams\ T[ N[ and Frederick\ W[ J[\ Kraft Recovery Boiler Physical and Chemical Processes[ America Paper Institute\ New York\ 0877[ 19[ Vanka\ S[ P[\ Block!implicit multigrid solution of NavierÐStokes

13[ 14[ 15[

16[

17[

18[

29[

20[

21[

22[

23[

24[

25[

26[ 27[ 28[

152

equations in primitive variables[ Journal of Computational Physics\ 0875\ 54\ 027Ð047[ Bohren\ C[ F[ and Hu}man\ D[ R[\ Absorption and Scatterin` of Li`ht by Small Particles[ John Wiley and Sons\ New York\ 0872[ Fiveland\ W[ A[\ Discrete!ordinates solutions of the radiative trans! port equation for rectangular enclosures[ ASME Journal of Heat Transfer\ 0873\ 095\ 588Ð695[ Carlso\ B[ G[ and Lathrop\ K[ D[\ Transport theory*the method of discrete ordinates[ Computin` Methods in Reactor Physics\ ed[ H[ Greenspan\ C[ N[ Kelber and D[ Okrent[ Gorden + Breach\ New York\ 0857[ Chandrasekhar\ S[\ Radiative Transfer[ Oxford at Clarendon Press\ London\ 0849[ Salcudean\ M[\ Nowak\ P[ and Abdullah\ Z[\ Cold ~ow com! putational model of a recovery boiler[ Journal of Pulp and Paper Science\ 0882\ 08"4#[ Matys\ P[\ Sabhapathy\ P[\ Nowak\ P[\ Abdullah\ Z[ and Salcudean\ M[\ E}ects of air and liquor delivery on carryover in a Kraft recovery furnace[ Paper presented at the 0884 International Chemi! cal Recovery Conference\ Toronto\ Ontario\ 13Ð14 April 0884[ Nowak\ P[\ Matys\ P[\ Sabhapathy\ P[\ Abdullah\ Z[ and Salcudean\ M[\ Numerical study of a Kraft recovery furnace[ Paper presented at the 0884 International Chemical Recovery Conference\ Toronto\ Ontario\ 13Ð14 April 0884[ Tse\ D[\ Gartshore\ I[ and Salcudean\ M[\ A numerical study of vorticity dynamics for a row of jets in a cross~ow[ Proceedin`s of the 04th Canadian Con`ress of Applied Mechanics "CANCAM 84#\ Vol[ 1\ pp[ 1\ 513Ð514[ Victoria\ British Columbia\ 17 MayÐ0 June 0884[ Abdullah\ Z[\ Salcudean\ M[\ Nowak\ P[\ Xiao\ Z[ Z[\ Savage\ M[\ Markovic\ C[\ Uloth\ V[ and Thom\ P[\ An initial validation of a cold ~ow mathematical model of a recovery boiler[ 0882 TAPPI Engineering Conference ^ published in TAPPI Journal\ May 0883\ 66"4#\ 038Ð046[ Quick\ J[\ Gartshore\ I[ and Salcudean\ M[\ Interaction of opposing jets with relevance to recovery furnaces[ Paper 5!3\ Ninth Sym! posium on Turbulent Shear Flows[ Kyoto\ Japan\ 05Ð07 August 0882[ Abdullah\ Z[\ Nowak\ P[\ Salcudean\ M[ and Gartshore\ I[ S[\ Investigation of interlaced and opposed jet arrangements for recov! ery furnaces[ Proceedin`s of the 0881 TAPPI En`ineerin` Confer! ence[ Boston\ Mass[\ 03Ð06 September 0881\ pp[ 092\ 011[ Quick\ J[\ Gartshore\ I[ S[ and Salcudean\ M[\ Interaction of oppos! ing jets with relevance to recovery furnaces[ Proceedin`s of the 0881 TAPPI En`ineerin` Conference[ Boston\ Mass[\ 03Ð06 September 0881\ pp[ 012Ð013[ Tse\ D[\ Gartshore\ I[ and Salcudean\ M[\ A study of multiple jets with attention to recovery boiler ~ow _elds[ Proceedin`s of the 0881 International Chemical Recovery Conference[ Seattle\ WA\ Book 0\ 6Ð00 June 0881\ pp[ 070Ð085[ Salcudean\ M[\ Nowak\ P[ and Abdullah\ Z[\ Mathematical mod! elling of recovery furnaces[ Proceedin`s of the 0881 International Chemical Recovery Conference[ Seattle\ WA\ Book 0\ 6Ð00 June 0881\ pp[ 086Ð197[ Gartshore\ I[ S[\ Aghdasi\ F[ and Salcudean\ M[\ Physical modelling of recovery boilers ] problems and prospects[ Paper presented at the 30st Canadian Chemical Engineering Conference\ Vancouver\ October 0880[ Salcudean\ M[\ Gartshore\ I[ S[\ Nowak\ P[\ Abdullah\ Z[\ Tse\ D[ and Quick\ J[\ Mathematical modelling of ~ow patterns in recovery boilers[ Paper presented at the 30st Canadian Chemical Engineering Conference\ Vancouver\ October 0880[ Gartshore\ I[ S[\ Salcudean\ M[ and Quick\ J[\ Bifurcation of ~ow pattern from opposing jets con_ned in a cavity[ Paper presented at the 02th CANCAM\ 0880[ Nowak\ Z[ P[ and Salcudean\ M[\ Turbulent ~ow calculations by the nonlinear multigrid method[ Akademie Verla`\ Z[ an`ew\ Math[ Mech[\ 0885\ 65"3#[ Tse\ D[\ Abdullah\ Z[\ Nowak\ P[\ Salcudean\ M[ and Gartshore\ I[ S[\ Evaluation of boundary conditions and calculations of jet ~ows in recovery furnaces[ Journal of Pulp and Paper Science\ 0885\ 11"0#[