Effect of Coil Design on the Temperature and Velocity Fields during Solidification in Electromagnetic Stirring Processes

Abstract

This paper examines the role of induction coil design on stirring of molten metal in electromagnetic (EM) solidification processes. A model is presented to describe the EM, heat transfer, and fluid flow phenomena in these processes. It is based on a dual-zone description of the mushy region, and accounts for damping of turbulence by the solidified crystallites. The electromagnetic field equations were solved using the mutual inductance technique, while the temperature and turbulent flow fields were calculated using the control volume method. Calculations were performed for solidification of an Al–Cu alloy placed in a stationary magnetic field generated by an induction coil. The effect of coil design on the flow structure was investigated for three different coil positions. It was found that changing the coil position significantly alters the flow pattern from four recirculating loops when the coil is above the midsection of the melt to two loops, typical of a travelling magnetic field, when the coil is at the base of the melt. This significantly modifies the rate of solidification across the ingot, as well as the temperature gradient, in the mushy region. The decay of the velocity and turbulent fields in the mushy region was found to be exponential, with the maximum rate of decay at the solidification front. These results indicate that through changes in coil design, it is possible to control the flow characteristics and solidification behavior in the molten metal.

1. Introduction

Electromagnetic (EM) stirring technologies have found extensive use in solidification processing of metals as a means to produce castings with a fine-grained equiaxed structure. The flow serves to minimize thermal gradients in the molten pool¹⁾ and rapidly dissipates the melt superheat,²⁾ which in turn promotes an early occurrence of the columnar-to-equiaxed transition.^3,4) In addition, it increases fragmentation of the secondary dendrite arms, leading to a higher nucleation potential in the melt.^1,5,6) Clearly, understanding the flow characteristics of an electromagnetically stirred melt in the mushy region is critical in order to control the grain structure of cast alloys.

Electromagnetic stirring is generally accomplished by induction from a time-varying magnetic field using a suitable induction coil. The passage of an alternating current induces eddy currents in the melt, which are dissipated as heat (Joule heating). Also, the interaction between the induced current and the applied magnetic field produces EM forces, which drive the flow in the melt. The flow behavior in EM stirred melts is generally controlled by the applied magnetic field used, namely stationary,^2,7,8) rotating,⁹⁾ or traveling magnetic fields.^5,10) Moreover, the resulting flow field is highly sensitive to variations in coil design, namely the coil configuration^11,12) and the location of the coil with respect to the melt.^13,14)

During past decades, investigations by Evans,^15,16) Fautrelle,^17,18) Szekely,¹⁹⁾ and El-Kaddah^10,20) led to a good understanding of the key flow phenomena in EM stirred melts. Their models involved the solution of the Maxwell and turbulent Navier-Stokes equations to describe the electromagnetic and velocity fields in the liquid metal. Implementation of this approach to solidification processes, however, is somewhat complicated due to the difficulties associated with the presence of the two-phase mushy region.

There are essentially two types of flow models for the mushy region, namely single-zone^7,21,22) and dual-zone models.^23,24) The single-zone model treats the entire two-phase region as a porous medium, with the flow being damped via Darcy’s law. The dual-zone model divides the mushy region into two domains: (i) the suspended particle region, and (ii) the fixed particle region, with the transition between regions occurring at the coherency point. In modeling EM solidification processes, the former model was used by Rappaz and coworkers⁷⁾ and Prescott and Incropera,²²⁾ while the latter was used by Pardeshi et al.²⁴⁾ In each of these studies, an approximate analytical solution was used to describe the EM force field in the system. This limits the capability of these models to address the role of the coil design on the solidification behavior.

Recently, the authors²⁵⁾ have proposed a new dual-zone model for EM solidification systems. This model includes numerical solution of the EM field, using the coil geometry and operating conditions as fundamental input parameters. The model also accounts for the damping of turbulence, a characteristic of EM stirred flows, in the two-phase mushy region due to interactions between the turbulent eddies and the solid crystallites. This paper investigates the influence of coil design on the flow characteristics for a single-phase induction coil.

2. Model Formulation

The mathematical description of the problem for any EM stirred solidification system is readily given in terms of the continuity, momentum, and heat transfer equations, along with appropriate subsidiary relationships defining the electromagnetic force field, Joule heating, and turbulence parameters in the system. This formulation will be given in broad outline only, since the derivation of the model equations are given in reference.²⁶⁾ Even further elaboration on the necessary steps in modeling such systems may be found in other references.^8,10,18)

2.1. Formulation of the Electromagnetic Field Problem

To allow for the presence of magnetic shields (chill block) and arbitrary coil configuration, it is convenient to use the mutual inductance method for two dimensional systems. In this method, the conducting region(s), with electrical conductivity σ, is divided into a number of elementary circuits with constant current density, J. The current density in each circuit is given by   

$$\begin{array}{l}\oint \mathbf{J}.d{\mathbf{l}}_i=\hfill \\ -j\omega \sigma \left[\sum _{m=1}^{metal}{M}_{i,m}{\left(\mathbf{J}\cdot \mathbf{S}\right)}_m+\sum _{c=1}^{shields}{M}_{i,c}{\left(\mathbf{J}\cdot \mathbf{S}\right)}_c+\sum _{k=1}^{coil}{M}_{i,k}{I}_k\right]\hfill \end{array}$$(1)

where S and l are the respective lengths and surface areas of the elementary circuits, j is the square root of –1, ω is the angular frequency, and M_i,j is the mutual inductance given by   

$$M_{i,k}=\frac{{\mu }_0}{4\pi }∯\frac{d{\mathbf{l}}_k\cdot d{\mathbf{l}}_i}{r\text{'}}$$(2)

From the current density, the magnetic flux density, B, Joule heating, q_em, and EM force field, F_em, are respectively determined by the following relations:   

$$\mathbf{B}=-\frac{\nabla \times \mathbf{J}}{j\omega \mathrm{\sigma }}$$(3)

  

$$S_{em}=\frac{1}{2}Re\left(\frac{\mathbf{J}\cdot {\mathbf{J}}^{*}}{\mathrm{\sigma }}\right)$$(4)

where the asterisk denotes the complex conjugate.

2.2. Formulation of the Heat Transfer Problem

Assuming that solid fraction, f_s, varies linearly with temperature, the differential equation representing the energy balance with phase change as a result of solidification may be written as:   

$$\rho C^{*}\frac{DT}{Dt}=\nabla \cdot \left(k^{*}\nabla T\right)+S_{em}$$(5)

where C^* is the effective specific heat given in terms of the specific heat c_p, latent heat L, and the solidus and liquidus temperatures, T_S and T_L, respectively, by   

$$C^{*}=\{\begin{array}{ccc}{c}_{p,l}& & T\ge T_L\\ c_{p,s}+\frac{L}{T_L-T_S}& & T_S<T<T_L\\ c_{p,s}& & T\le T_L\end{array}$$(6)

In Eq. (6), the l and s subscripts correspond to values for the liquid and solid phases, respectively. In convective turbulent flows, the effective thermal conductivity, k^*, in the bulk liquid, mushy region, and solidified metal is given by:   

$$k^{*}=\{\begin{array}{l}\\ \\ \\ \end{array} \begin{array}{ccc}{k}_l+\frac{{\mu }_t{c}_{p,l}}{Pr}& & T\ge T_L\\ \left(1-f_S\right)k_l+f_S{k}_s& & T_S\le T<T_L\\ k_s& & T\le T_S\end{array}$$(7)

For bulk liquids (T>T_L), the effective thermal conductivity is given by the sum of the turbulent and molecular thermal conductivities, where μ_t and Pr are the turbulent viscosity and turbulent Prandtl number, respectively. In the mushy region, the conductivity is a weighted sum of the solid and effective liquid conductivities.

2.3. Formulation of the Fluid Flow Problem

For the three flow domains, namely the bulk liquid, suspended particle, and fixed particle regions, the continuity and Navier-Stokes equations can be generally expressed as:

Continuity:   

$$\nabla \cdot \mathbf{u}=\text{0}$$(8)

where   

$$\mathbf{u}=f_s{\mathbf{u}}_s+\left(1-f_s\right){\mathbf{u}}_l$$(9)

Navier-Stokes:   

$$\rho \frac{\partial \mathbf{u}}{\partial t}+\rho \mathbf{u}\cdot \nabla \mathbf{u}=-\nabla P+\mu {\nabla }^2\mathbf{u}+\lambda \left({\mathbf{F}}_{em}+{\mathbf{F}}_d\right)-\frac{{\mu }_l}{K}\left(\mathbf{u}-{\mathbf{u}}_s\right)$$(10)

where F_em, F_d, K, and λ represent the bulk electromagnetic force, turbulent damping force, permeability, and a switch parameter, respectively. In the bulk liquid and suspended particle regions, λ is equal to unity and is zero in the fixed particle region. Since it is assumed that the velocity of solid is equal to velocity of liquid in the bulk and suspended particle region, the permeability term is finite only in the fixed particle region. In this study, the permeability was determined using the Carman-Kozeny equation.²⁷⁾

In this work, it is assumed that the flow in the fixed particle region is laminar. This assumption is reasonable, as the interparticle spacing in this region is on the order of the Kolomogrov turbulence length scale. For the bulk liquid and suspended particle regions, the turbulent characteristics of the flow are defined by the low-Re k-ε model.²⁸⁾ The governing equations for the turbulent kinetic energy, k, and the turbulent energy dissipation, ε, are:

k-equation   

$$\rho \frac{\partial k}{\partial t}+\rho \mathbf{u}\cdot \nabla k=\nabla \cdot \left(\frac{{\mu }_t}{P_k}\nabla k\right)+G_{ij}-\rho \mathit{\epsilon }$$(11)

ε-equation   

$$\rho \frac{\partial \mathit{\epsilon }}{\partial t}+\rho \mathbf{u}\cdot \nabla \mathit{\epsilon }=\nabla \cdot \left(\frac{{\mu }_t}{P_{\mathit{\epsilon }}}\nabla \mathit{\epsilon }\right)+C_1{ }^{*}{G}_{ij}\frac{\mathit{\epsilon }}{k}-C_2{ }^{*}\rho \frac{{\mathit{\epsilon }}^2}{k}$$(12)

where C₁^* and C₂^* are constants. The turbulent viscosity is defined by   

$${\mu }_t=C_{\mu }{ }^{*}\rho \frac{k^2}{\mathit{\epsilon }}$$(13)

where C_μ^* is a constant. The effective viscosity is the sum of the laminar and turbulent components. The rheological viscosity in the suspended particle region is represented using the correlation developed by Thomas.²⁹⁾

2.4. Turbulent Damping Model

The turbulent damping force in the suspended particle region is   

$${\mathbf{F}}_d=\left(1-f_L{ }^2\right)\nabla .{\overline{\tau }}_t$$(14)

where f_L is the liquid fraction and τ_t is the turbulent Reynolds stress. It is assumed that the turbulent stress is isotropic and proportional to the turbulent kinetic energy. Therefore, the divergence term in Eq. (14) is equal to   

$$\nabla .{\overline{\tau }}_t=c_d{ }^{1/2}\rho \nabla k$$(15)

where c_d is a constant equal to 0.08.³⁰⁾ In Eq. (15), the turbulent grad term is equal to the difference in turbulent kinetic energy between the liquid and the crystallite surface, where k is taken to be zero. For a simple cubic arrangement of crystallites of size d with distance between vertices D_f, the grad term may be approximated by:   

$$\nabla k=\frac{2k}{\left(D_f-d\right)}{\stackrel{⌢}{\mathbf{e}}}_j$$(16)

To determine the instantaneous particle size, the solid fraction of any particular cell in the suspended particle region is given by   

$$f_s=\frac{\pi d^3\beta }{6D_f^3}$$(17)

where β is the sphericity, defined as the ratio of the particle volume to the volume of the spherical envelope of the crystallite. For spherical particles, β is equal to unity. From Eq. (17), the instantaneous diameter can be written as:   

$$d=\gamma D_f{f}^{1/3}\hspace{1em}\mathrm{\hspace{0.17em} }\text{;}\hspace{0.17em}\mathrm{\hspace{0.17em} }\hspace{1em}\gamma ={\left(\frac{6}{\pi \beta }\right)}^{1/3}$$(18)

From Eqs. (14), (15), (16) and (18), the damping force can be given in terms of solid fraction by:   

$${\mathbf{F}}_d=\frac{2\left(f_s{ }^2-2f_s\right)\sqrt{c_d}\rho k}{D_f\left(1-\gamma f_s{ }^{1/3}\right)}{\hat{\mathbf{e}}}_j$$(19)

It should be noted that in previous formulations²⁵⁾ the turbulent damping force was given in terms of the coherency solid fraction. In this formulation, the coherency is tied to the sphericity such that the coherency is equal to πβ/6, and for spherical particles is equal to 0.521.

3. Solution Technique

To calculate the electromagnetic field, the metal and adjacent conducting regions were divided into a number of elementary circuits. A finer grid was used near the outer radius of the conducting domains to accurately describe the decay of the magnetic field in the skin depth. The current density in each elementary circuit was calculating solving the set of J equations using Choleski factorization method.³¹⁾ Upon obtaining the J field, the magnetic flux density B was evaluated from Ampere’s law.³²⁾ These two fields were then used to calculate the electromagnetic force field and Joule heating.

The control volume technique³³⁾ was used to solve the equations for the temperature, fluid flow, and turbulent fields. In this technique, the solution domain was divided into a number of contiguous volume elements. The equations were discretized over each individual control volume, and the resulting discretized equations were solved iteratively using tri-diagonal algorithm technique (TDMA). In the calculation of the velocity field the SIMPLE algorithm was used.³⁴⁾

4. Results and Discussion

In this section, we shall present the effect of coil design on the computed velocity and temperature fields during solidification. This study focuses on the influence of the coil position on the flow and solidification phenomena. Calculations were carried out for unidirectional solidification of Al–Cu alloy in a bottom chill mold surrounded by an induction coil, as shown in Fig. 1. The solidus and liquidus temperatures of the alloy are 564 and 648°C, respectively. The inner diameter of the mold and the ingot height are both 80 mm, respectively, and the thickness of the stainless steel chill block is 80 mm. The coil used in these calculations consists of three turns, each with a diameter of 130 mm and a spacing of 10 mm between each turn. Three coil positions were examined, and were chosen such that the coil center was aligned with the top, middle, and bottom of the melt. Stirring of the molten metal during solidification was carried out for a coil current of 200 A, and a frequency of 500 Hz. In this study, the crystallites were assumed to be spherical (β=1), which corresponds to globular solidification morphology. The final grain size, D_f, was taken to be 200 μm.

Fig. 1.

Sketch of unidirectional EM solidification system.

Figures 2, 3, 4 show the initial flow fields for the three coil positions investigated. These figures show that the coil location has a substantial effect on the flow pattern. The flow corresponding to the middle coil position is characterized by two recirculating loops of equal size, Fig. 2. Figures 3 and 4 show that the flow pattern changes drastically when the coil is placed at either the top or bottom of the melt, respectively. When the coil is located at the top, the lower loop grows considerably in size at the expense of the upper loop, while the flow for the bottom coil position is characterized by essentially one loop similar to those typically found for traveling magnetic fields.¹⁰⁾ Here it should be mentioned that when the coil is above the midplane of the melt, the flow in the upper loop is clockwise, and is counterclockwise in the lower loop. For the bottom coil position, the flow is opposite that of the middle and top positions at the metal/chill interface. These results suggest that the coil location can be used to control the flow directionality in the melt.

Fig. 2.

Initial flow pattern for the middle coil location.

Fig. 3.

Initial flow pattern for the top coil location.

Fig. 4.

Initial flow pattern for the bottom coil location.

The coil position affects the turbulent characteristics of the flow as well. Figures 5, 6, 7 show the turbulent viscosity, normalized with respect to the laminar viscosity (μ^*=μ_t /μ_l), for the middle, top, and bottom coil positions, respectively. The turbulent field for the middle coil position, Fig. 5, is essentially symmetrical, with the maximum values of turbulent viscosity being found at the vortex centers and the free surface, where the shear strain rates are highest. By placing the coil at the top of the metal, where the flow is characterized by two unequal recirculating loops, the turbulent viscosity is higher than that for the middle coil position. Furthermore, the maximum values of turbulent viscosity are located in the vortices of each loop, Fig. 6. When the coil is placed at the base of the melt, where the flow is characterized by a single loop, the spatial variation of the turbulent viscosity is completely different compared to that of the other two coil positions. Turbulence is stronger than that of the other cases, and is more or less uniform throughout the melt, Fig. 7. Since the flow affects the temperature field in both the liquid and mushy regions, these results suggest that the coil location has a strong influence on the solidification rate.

Fig. 5.

Initial turbulent field for the middle coil position.

Fig. 6.

Initial turbulent field for the top coil position.

Fig. 7.

Initial turbulent field for the bottom coil position.

The changes in solidification rate with respect to the coil position are illustrated in Figs. 8, 9, 10. These figures show the velocity and temperature fields after 120 s of solidification for the middle, top, and bottom coil positions, respectively. It is seen that stirring changes the solidification behavior from unidirectional to non-uniform, as indicated by the radial variation of temperature in the mushy region. In addition, the local rate of solidification depends on the direction of the flow. Inspection of these figures shows that when the flow travels away from the chill plate, the rate of solidification is higher, and is depressed when the flow is travelling downward as a result of remelting the solidification front.

Fig. 8.

Computed velocity and temperature fields in the suspended particle region for the middle coil position after 120 s.

Fig. 9.

Computed velocity and temperature fields in the suspended particle region for the top coil position after 120 s.

Fig. 10.

Computed velocity and temperature fields in the suspended particle region for the bottom coil position after 120 s.

The dependence of solidification rate on the flow field is reflected in the obtained changes of solidification rate across the melt with coil location. For the top coil position, the maximum rate of solidification is near the mold wall, while it is at the centerline of the metal for the bottom coil position. This point is further illustrated in Fig. 11, which shows the radial temperature variation for the top and bottom coil positions at z=20 mm. As seen in this figure, the temperature at the centerline is at a maximum for the bottom coil position and decreases across the melt. For the top coil position, the temperature is minimum at the centerline and increases along the radius.

Fig. 11.

A comparison of the radial temperature distributions for the top and bottom coil positions after 120 s (z=20 mm).

Since the convective heat flux depends on the velocity and turbulent intensity, it is instructive to examine the effect of coil location on the decay of the flow field. As seen in Figs. 8, 9, 10, the coil position has little effect on the magnitude of the velocity in the suspended particle region. The decay in the magnitude of the velocity during solidification is shown in Fig. 12. As seen in this figure, the flow in the suspended particle region decays rapidly in the initial stage of solidification, which corresponds to solid fraction less than 0.1, and the rate of decay decreases as solidification progresses. Furthermore, the rate of decay is essentially independent of the coil position, except near the interface between the suspended and fixed particle regions.

Fig. 12.

Decay of velocity in the suspended particle region for top and bottom coil positions (r=20 mm, z=30 mm).

Figure 13 shows the corresponding decay of the turbulent viscosity normalized with respect to the laminar value. This figure shows that the decay of the turbulent viscosity is similar to that of the velocity shown in Fig. 12. From these results, it is apparent that the coil position has little or no effect on the overall rate of solidification in the metal. These findings are consistent with the cooling rate measurements by Campanella et al. for different coil positions.⁶⁾

Fig. 13.

Decay of turbulent-to-laminar viscosity ratio in the suspended particle region for top and bottom coil positions (r=20 mm, z=30 mm).

5. Conclusions

A mathematical model for EM stirred solidification processes has been developed. The presented solidification model is based on the numerical solution of the electromagnetic field equations using the mutual inductance method, and extends the dual-zone velocity and temperature field models to incorporate turbulent flow. The model also accounts for the damping of turbulence at the crystallite surfaces. The model provides a vivid description of the evolution of the various electromagnetic, heat transfer, and flow phenomena in the system during solidification. Furthermore, the model contains the flexibility such that it can describe these phenomena for any coil configuration, position, and operating conditions, making it a powerful tool in the analysis and/or design of EM solidification processes.

The model was used to investigate the effect of coil location on the flow and temperature evolution during solidification in a stationary magnetic field. The results of this study shows that the coil position significantly modifies the flow pattern in the melt, and alters the directionality of the flow in the bulk liquid. The coil location was also found to affect the local solidification rate. It was also found that the decay of the velocity and turbulent viscosity primarily occurs in the earliest stages of solidification (f_s < 0.1), and is independent of coil position. It is hoped that the model presented in this work will provide a sound theoretical framework and will allow for a more in-depth understanding of EM stirring during solidification and for precise control of the flow structure in EM stirred systems.

References

• 1)   B.  Willers,  S.  Eckert,  U.  Michel,  I.  Haase and  G.  Zouhar: Mater. Sci. Eng. A, 402 (2005), 55.
• 2)   C.  Vivés and  C.  Perry: Int. J. Heat Mass Transf., 29 (1986), 21.
• 3)   B.  Chalmers: J. Aust. Inst. Met., 8 (1963), 255.
• 4)   J. A.  Dantzig and  M.  Rappaz: Solidification, EPFL Press, Lausanne, (2009), 458.
• 5)   W.  Griffiths and  D.  McCartney: Mater. Sci. Eng. A, 216 (1996), 47.
• 6)   T.  Campanella,  C.  Charbon and  M.  Rappaz: Metall. Mater. Trans. A, 35 (2004), 3201.
• 7)   T.  Campanella,  C.  Charbon and  M.  Rappaz: Scr. Mater., 49 (2003), 1029.
• 8)   N.  El-Kaddah and  J.  Szekely: J. Fluid Mech., 133 (1983), 37.
• 9)   V.  Travnikov,  P. A.  Nikrityuk,  K.  Eckert,  D.  Räbiger,  S.  Odenbach and  S.  Eckert: IOP Conf. Ser.: Mater Sci. Eng., 33 (2012), 012048.
• 10)   N.  El-Kaddah,  J.  Szekely and  G.  Carlsson: Metall. Trans. B, 15 (1984), 633.
• 11)   J. W.  McKelliget and  N.  El-Kaddah: J. Appl. Phys., 64 (1986), 2948.
• 12)   V.  Bojarevics,  K.  Pericleous and  M.  Cross: Metall. Mater. Trans. B, 31 (2000), 179.
• 13)   V.  Bojarevics,  A.  Roy and  K. A.  Pericleous: Magnetohydrodyn., 46 (2010), 317.
• 14)   C.  Stelian,  Y.  Delannoy,  Y.  Fautrelle and  T.  Duffar: J. Cryst. Growth, 275 (2005), 1571.
• 15)   E. D.  Tarapore and  J. W.  Evans: Metall. Trans. B, 7 (1976), 343.
• 16)   E. D.  Tarapore,  J. W.  Evans and  J.  Langfeldt: Metall. Trans. B, 8 (1977), 179.
• 17)   Y.  Fautrelle: J. Fluid Mech., 102 (1981), 405.
• 18)   N.  El-Kaddah,  J.  Szekely,  E.  Taberlet and  Y.  Fautrelle: Metall. Trans. B, 17 (1986), 687.
• 19)   K.  Nakanishi,  T.  Fujii and  J.  Szekely: Ironmaking Steelmaking, 2 (1975), 193.
• 20)   J. L.  Meyer,  J.  Szekely and  N.  El-Kaddah: Trans. Iron Steel Inst. Jpn., 27 (1987), 25.
• 21)   R.  Mehrabian,  M.  Keane and  M. C.  Flemings: Metall. Trans., 1 (1970), 455.
• 22)   P. J.  Prescott and  F. P.  Incropera: J. Heat Transf., 117 (1995), 716.
• 23)   C. M.  Oldenburg and  F. J.  Spera: Numer. Heat Transf. B, 21 (1992), 217.
• 24)   R.  Pardeshi,  A. K.  Singh and  P.  Dutta: Numer. Heat Transf. A, 55 (2008), 42.
• 25)   G. M.  Poole and  N.  El-Kaddah: J. Iron Steel Res., Int., 19, Suppl. 1, (2012), 252.
• 26)   G. M.  Poole and  N.  El-Kaddah: Metall. Trans. B, 44 (2013), No. 6, 1531.
• 27)   P. C.  Carman: Trans. Inst. Chem. Eng., 15 (1937), 150.
• 28)   W. P.  Jones and  B. E.  Launder: Int. J. Heat Mass Transf., 15 (1972), 301.
• 29)   D. G.  Thomas: J. Colloid Sci., 20 (1965), 267.
• 30)   B. E.  Launder and  D. B.  Spalding: Mathematical Models of Turbulence, Academic Press, London, (1972), 78.
• 31)   N. V.  Kopchenova and  I. A.  Maron: Computational Mathematics, Mir Publishers, Moscow, (1975), 63.
• 32)   J. D.  Jackson: Classical Electrodynamics, 2nd ed., Wiley, New York, (1975), 173.
• 33)   W. M.  Pun and  D. B.  Spalding: Report HTS/76/2, Imperial College, London, (1972).
• 34)   S. V.  Patankar: Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing, New York, (1980).