ASPECTS OF MOULD DESIGN IN ELECTROSLAG CASTING by JAYANT MORESHWAR SATHAYE B.Tech., Indian Institute Of Technology, Madras, 1977 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department Of Metallurgical Engineering We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA February 1983 © Jayant Moreshwar Sathaye, 1983 In presenting t h i s thesis in p a r t i a l fulfilment of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t freely available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of Metallurgical Engineering The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Abstract Electroslag casting i s a process technique used to produce castings using electroslag melting. It i s l i k e l y that channel cooled Aluminium moulds w i l l be employed. Thus, the c r i t e r i a for the design of moulds to be used in t h i s process were examined. Temperatures in a mould section were measured and compared to those calculated from a model developed assuming 3-D quasi steady state heat flow and reasonable agreement was found. In the model the two major boundary conditions are the heat transfer c o e f f i c i e n t to the coolant and the heat flux on the hot face. The former can be estimated from established correlations for heat transfer in tubes such as the Petukhov and the Sieder-Tate r e l a t i o n s . The l a t t e r had to be measured and a heat flux sensor was developed for t h i s . Proposals were made to q u a l i t a t i v e l y explain the heat flux v a r i a t i o n along the mould height, but i t was not possible to estimate the heat flux for a casting of a r b i t r a r y shape in a simple manner. According to the slag c r y s t a l l i s a t i o n mechanism proposed the influence of mould on the surface qua l i t y of the casting was found to be n e g l i g i b l e . Table of Contents Abstract i i L i s t of Tables v L i s t of Figures v i Acknowledgements ix Chapter I INTRODUCTION 1 1.1 The Electroslag Process 1 1 .2 The Mould 2 1.2.1 Heat Flow In And To The Mould 2 1.2.2 Heat Transfer To The Mould In Mathematical Models 7 1.2.3 Design Of Moulds 8 1.2.4 Water Cooled Moulds Used Elsewhere 11 1.2.5 Probable Mould Type In ESC 11 1.3 Surface Quality 12 1.3.1 Slag Skin C r y s t a l l i s a t i o n 13 1.3.2 Slag Skin E f f e c t s On Surface Quality 15 1 . 4 Objectives 16 Chapter II EXPERIMENTAL 17 2.1 Experiments On The Furnace 17 2.1.1 The Electroslag Furnace At UBC 17 2.1.2 Temperature Measurements 18 2.1.3 The Heat Flow Sensor 22 2.1.4 V e r i f i c a t i o n Of Sensor Data 22 2.1.5 Experimental Results 24 2.2 Observations On Slag Skins 36 2.2.1 Microstructure 36 2.2.2 Observations On The SEM 42 2.2.3 The 70:30 Slag Cap Near The Liquid Metal Meniscus 48 Chapter III CALCULATIONS 50 3.1 Calculations In 2 Dimensions .50 3.2 Calculations In 3 Dimensions 52 3.3 Speculation On The Effe c t Of Grooves On Surface Quality 60 Chapter IV DISCUSSION 62 4.1 Failur e Of 2-D Model 62 4.2 Apparent Success Of 3-D Quasi Steady State Model ..65 4.3 Factors Af f e c t i n g The Heat Flux To The Mould 67 4.4 C r y s t a l l i s a t i o n In Slag Skins 72 4.5 Slag Composition E f f e c t s On The Heat Flux 74 4.6 Factors Af f e c t i n g Surface Quality 76 iv Chapter V SUMMARY 78 REFERENCES 79 APPENDIX A - SLAG NOTATION 83 APPENDIX B - MISC. DATA AND CORRELATIONS 84 APPENDIX C - FORTRAN PROGRAM TO CALCULATE TEMPERATURES IN CHANNEL COOLED MOULDS 86 APPENDIX D - FORTRAN PROGRAM REFERRED TO IN SECTION 3.3 .101 APPENDIX E - F.E. FORMULATION FOR STEADY HEAT TRANSFER IN A PLANE' 117 V L i s t of Tables I. Thermal Resistance - from Kusamichi, et a l .....6 I I . Record of experimental runs 25 I I I . Comparison: Heat flow and Power input 35 IV. Chemical analysis of some skins 48 V. Possible Error in heat flow measurement 65 VI. Temp, drop through slag skin 67 VII. Heat flow through a i r gap 68 v i L i s t of Figures 1. Heat Loss vs. Parameter F 4 2. Heat Flux Variation 4 3. Sketch of Mould Assembly 19 4. Views of Mould Assembly 19 5. Dimensions of Mould Section used 20 6. Dimensions of Copper Section 21 7. Sketch of Heat Flux Sensor 23 8. For V e r i f i c a t i o n of Sensor 23 9. Electrode used in Run 10 26 10. Temp. History of RUN 3 • 26 11. Smoothened Temp. History of Run 3 27 12. Calculated Temp, for Run 3 27 13. Temp.History of Run 2 28 14. Calculated Temp, for Run 2 28 15. Temp.History of Run 6 29 16. Calculated Temp, for Run 6 29 17. Temp.History of Run 10 30 18. Calculated Temp, for Run 10 30 19. Temp.History of Run 8 31 20. Heat Flux Measurement-Run 2 31 21. Heat Flux Measurement-Run 3 32 22. Heat Flux Measurement-Run 4 ....32 23. Heat Flux Measurement-Run 6 33 24. Heat Flux Measurement-Run 11 33 25. Sample comparison between Calculated temperature and v i i measured 34 26. SEM picture of 70:30 skin ...37 27. Spike on 70:30 skin 37 28. 70:15:15 slag skin 38 29. 50:30:20 slag skin 38 30. Skin of S i l i c a , Magnesia containing slag 39 31. Skins of complex slag 39 32. Industrial 70:30 slag skin 40 33. Industrial 50:30:20 slag skin 40 34. 70:30 slag cap near metal meniscus 41 35. 70:30 slag cap near top 41 36. EDX analysis of 70:30 skin 44 37. EDX analysis of 70:15:15 skin 45 38. EDX analysis of 50:30:20 skin 46 39. EDX analysis of complex slag skin 47 40. Temp, f i e l d c a l c . with 2-D assumptions 51 41. The s l i c e which was d i s c r e t i s e d 51 42. Calculated Temp, f i e l d s 56 43. Calculated Temp, f i e l d s 57 44. Calculated Temp, f i e l d s 58 45. Section taken for F.E. cal c u l a t i o n 58 46. Calculated results of speculation 59 47. Error in heat flux measurement 64 48. Hypothetical heat flux v a r i a t i o n 64 49. Phase diagram for CaF 2-AI2O3 69 50. Phase diagram for CaF 2-Al 2 03-CaO 73 51. Shape functions-Rect. isoparametric element 118 v i i i 52. Shape functions-Triangular element i x Acknowledgement I would l i k e to thank Prof. Alec M i t c h e l l for a l l the time he spared for me and e s p e c i a l l y for h i s s t y l e of guidance. Gus S i d l a ' s assistance was invaluable and i s very much appreciated. Ed Barry and Rudy Cardeno were very h e l p f u l during the experimental stages. The pleasant nature of the s t a f f in the department and of fellow graduate students made my stay very enjoyable. I also wish to express my gratitude to Consarc Corporation, Rancocas, New Jersey for t h e i r f i n a n c i a l support. 1 I. INTRODUCTION 1.1 The Electroslaq Process The u t i l i s a t i o n of the electroslag process began in the early 1940's with the Kellogg e l e c t r i c ingot process. The process was o r i g i n a l l y developed by R.K.Hopkins, but western industry did not involve i t s e l f seriously t i l l the 1960's. Investigators in the Soviet Union, meanwhile, discovered that the weld metal in an electroslag weld was purer and had fewer inclusions than the base metal. They then developed the process with the prime aim of improving the quality of air-melted e l e c t r i c furnace s t e e l . Since then, i t has been used to produce castings too. Following t h i s lead the process has been used in many other countries in the manufacture of ingots. A Japanese firm uses t h i s process to manufacture reformer tubes but apart from th i s there are almost no s i g n i f i c a n t engineering applications of electroslag cast metal, as i s , in the western world. The e l e c t r o s l a g process involves the melting of a s o l i d electrode through a l i q u i d e l e c t r i c a l l y conducting slag into a water cooled metallic mould. An e l e c t r i c potential i s applied across the electrode and the mould; The slag forms the only s i g n i f i c a n t resistance and almost a l l the heat i s generated within i t . The electrode i s melted progressively and the ingot s o l i d i f i e s in similar fashion, i . e . progressively and to a 2 certain extent d i r e c t i o n a l l y . It reduces, s i g n i f i c a n t l y , the number of inclusions and can control the s o l i d i f i c a t i o n rate thus eliminating shrinkage porosity and minimising segregation. The produced metal i s considered more r e l i a b l e . Attempts are being made to est a b l i s h e l e c t r o s l a g castings in engineering use. The term 'casting' implies that complex shapes, r e l a t i v e to ingots, w i l l be considered. This s i g n i f i e s that mould design w i l l be an important aspect in deliberations about electroslag casting and fresh designs w i l l be required for every new shape. 1.2 The Mould The mould, in general, i s a container for the l i q u i d metal and serves to shape the casting. Additionally i t should minimise the p r o b a b i l i t y of such defects as hot tears. Basic p r i n c i p l e s of mould design u t i l i s e d in conventional casting are applicable here. The extra feature i s that the mould must be able to abstract the heat fluxes l i a b l e to occur in the elec t r o s l a g process and should foster the formation of a smooth surface. 1.2.1 Heat Flow In And To The Mould Most published work on heat transfer to the mould in the electroslag process has concentrated on the conditions occurring in ingot remelting. Nearly 35-50 percent of the heat generated in the slag passes d i r e c t l y to the mould and another 20-25 3 percent passes i n d i r e c t l y to the mould. 1' 2' 3'" The heat flux has been found to r e f l e c t the pattern of heat generation with the peak heat flux being in the range of 1.6 MW/m2. Early work on the heat balance of the process was done by Soviet workers at Kiev, the essence of which i s reviewed in the book edited by Medovar et a l . 1 Depending on electrode configuration ( b i f i l a r , 3-phase, e t c . ) , the slag depth and the metal slag i n t e r f a c i a l area; between 46 to 65% of the power input goes d i r e c t l y to the slag. An empirical r e l a t i o n has been developed for t h i s as shown in f i g . 1. The heat flux too varies with the above. In addition, i t varies with the height of the ingot casting, the power input and the mould shape. Many measurements of heat flux are reported. In some cases two peaks appear, one at the a i r slag interface and the other at the slag metal interface and in other cases a peak does not occur, but there i s a broad band of high heat flux from the l i q u i d slag, as i s shown in f i g 2. The heat flux from the slag reduces with increasing l i q u i d metal head, increasing distance from the electrode and decreasing power input. As the slag skin thickness increases, the mid slag region heat f l u x / metal head flux r a t i o tends to come closer to 1 . Much of th i s has been confirmed by Holzgruber 3 and by J o s h i . 5 Similar energy balances and heat flows can be deduced from the work described below. Pocklington and P a t r i c k 6 report on heat flow in 3 phase slab production. The moulds were jacketed and 1 D approximation was used to calculate heat fluxes from thermocouples inserted in 90 80 H 40 • X / » 20-1 2 0 3 0 4 0 5 0 F = s l ° S - m o u l d .n ler foce a r e a s l a g - m e t a l interface a r e a " F i g u r e 1 - H e a t L o s s v s . P a r a m e t e r F i g u r e 2 - He a t F l u x V a r i a t i o n 5 the wall. Heat flux variations were found to depend on slag composition and the remelting conditions. Peaks of temperature were found in what they c a l l e d 'Form' 4 curves, but t h i s did not translate into peaks of heat fl u x . Maximum heat fluxes were of the order of 1 - 1.2MW/m2 for most slags and about 1.8 MW/m2 for 70:15:15 slag. Greater va r i a t i o n of heat flow around the periphery was found compared to conticast moulds. It was suggested that small heat flux d i f f e r e n t i a l s may minimise d i s t o r t i o n of the mould. Some experiments were c a r r i e d out by the Albany Metallurgy Research Centre. 7 A 3 piece jacketed mould was used in a laboratory furnace. They observed 3 peaks in the heat flux curve with maximum flux of about 3 MW/m2. The temperature d i f f e r e n t i a l of the cooling water between the i n l e t and the outlet was used to determine the heat flux and i t i s l i k e l y that very high precision measurements were required. Kusamichi et a l 8 conducted interesting experiments in which slag was frozen onto a water cooled copper pipe in a crucible of molten slag. Two thermocouples were kept very close to the pipe so that the slag s o l i d i f i e d around them and the measured temperatures were used to estimate the ef f e c t i v e thermal conductivities of slag skins, the thermal resistances, and the heat transfer c o e f f i c i e n t s . The resistance of the slag skin and that of the slag skin - mould interface were found to be of major importance and this i s shown in Table I. Differences were noted for d i f f e r e n t slags, but i t must be pointed out that c r y s t a l l i s a t i o n i s d i f f e r e n t in the real process. 6 Table I - Thermal Resistance - from Kusamichi et a l a. in Crucible Slag Thermal Resistance (s.deg./Cal) CaF 2 A1 20 3 Cao Total Bath- Crust Crust- Cu Pipe-Crust Pipe Pipe water 0 50 50 35.4 4.6 28.5 0.7 0.01 2.2 50 25 25 30.3 3.6 14.3 10.9 0.01 1.7 100 0 0 21 .7 2.9 5.4 12.0 0.01 1 .3 b. in ESR 0 50 50 0.81 0.19 0.48 0.09 0.007 0.04 50 25 25 0.43 0.09 0.21 0.10 0.007 0.02 In (b) read mould for pipe. Kondo, Kodama et a l , 9 measured heat fluxes of up to 0.8 - 1 MW/m2 in 800 mm i n d u s t r i a l moulds. A peak in the heat flux d i s t r i b u t i o n was seen at the slag metal interface l e v e l for 70:15:15 slag and i t was claimed that remelting of the slag skin of a 70:30 slag was observed. The difference between Soviet investigators and the others i s noteworthy. The Soviets developed a dire c t heat sensor which was thermally insulated from the rest of the mould, the signal being generated by a thermopile. The others merely inserted thermocouples into the mould wall and^ j u s t i f i e d the technique by suggesting that 1-D approximation was acceptable. 7 1.2.2 Heat Transfer To The Mould In Mathematical Models It i s of interest to examine the boundary conditions at the slag mould interface employed by the various mathematical models around. Many models did not include the slag, but one may compare the conditions at the top of the ingot. Sun and Pridgeon 1 0 determined the conductance, h, by plunging Copper rods into molten slag and metal and measuring heating rates. Their value for "h" to the mould varied from .006 to .00745 cals./cm 2sK. Carvajal and Geiger 1 1 used 0.4 cals./cm 2sK and they did not model the slag. E l l i o t t and Maulvault 1 2 used values from .01 to .02 cals./cm 2sK as did Paton and co-workers. 1 3 J o s h i 5 experimentally estimated i t to be .0115 cals./cm 2sK while B a l l a n t y n e 1 0 used .02. More advanced models including f l u i d flow have been developed by Kreyenberg and Schwerdtfeger 1 5 and by Dilawari, Choudhary and S z e k e l y . 1 7 ' 1 8 ' 1 9 These workers did not consider heat transfer c o e f f i c i e n t s but fixed the slag boundary to be at i t s liquidus temperature. In a recent paper Choudhary and S z e k e l y 1 9 used wall functions to estimate the heat transfer c o e f f i c i e n t to the slag skin boundary and predicted that the heat transfer c o e f f i c i e n t w i l l remain near constant. This s i t u a t i o n i s attributed to a nearly invariant temperature gradient close to the mould wall in the slag along the length of the mould face. In general, these modellers were interested in the pool shape and the l o c a l s o l i d i f i c a t i o n time at the centre of the ingot. 8 1.2.3 Design Of Moulds In the case of ingot remelting most moulds are jacketed and water cooled as a non-boiling system using p r i n c i p l e s as detailed by M i t c h e l l and Smailer. 2 Channel-cooled moulds (moulds with holes d r i l l e d in the walls) have also been used e s p e c i a l l y in the moving-mould and ingot-withdrawal techniques. Early moulds were made of steel( with spray cooling) but practice soon switched to copper. Steel moulds had to be r e l a t i v e l y thinner for good a x i a l heat conduction and the danger of burn-through was heightened. Further, copper moulds were found to l a s t much longer and so were cheaper in the long run. Larger moulds are made from slabs of thicker section since the stresses the mould must take are much higher. Cremisio and Z a k 2 0 ' 2 1 have discussed the design of slab moulds in d e t a i l . Mould design starts with the minimum water ve l o c i t y required to maintain non-boiling conditions. The mould wall thickness i s generally selected from experience and the rest of the detailed design follows from t h i s . They are of the view that f l a t surfaces in moulds face more demanding conditions than curved surfaces as in round moulds. The design should permit movement in the horizontal and v e r t i c a l planes while maintaining r e s t r a i n t in the t h i r d plane. Creep resulting from thermal stresses i s the chief problem in such moulds. Channel cooled moulds have s u f f i c i e n t r e s t r a i n t s b u i l t in as a result of their physical construction. However, they face higher temperatures on the hot face. Such a design may frequently not be used in f u l l - l e n g t h moulds since t h i s may 9 e n t a i l inordinately high water pressures to ensure s u f f i c i e n t l y high water-flow rates. Soviet researchers were the early developers of the process and have done extensive work in the f i e l d . The g i s t of t h i s work has been summarised in the book edited by Medovar et a l . 1 Highly involved empirical design techniques have been developed and i t appears that their work i s the most detailed published so far. Various methods of water cooling were investigated, e.g., jacket cooling, channel cooling and " r i b " cooling. Schemes for cooling with l i q u i d metals and organic l i q u i d s were formulated but i t i s not known i f t h i s was put into practice. Both evaporative and forced convection were considered with the l a t t e r being generally preferred. The Soviet approach to design i s from two d i r e c t i o n s ; from the slag to the mould and from the coolant to the mould. Using the empirical methods mentioned in section 1 . 2 . 1 the heat flux to mould wall i s estimated. Taking a safety factor the water ve l o c i t y i s chosen depending on the mould type ( i . e . the " e f f i c i e n c y " ) and from t h i s the detailed hydraulic design i s carried out. The expected maximum mould temperature at the mould-coolant interface i s estimated from the heat flux to the coolant, the coolant temperature and the heat transfer c o e f f i c i e n t . The heat transfer c o e f f i c i e n t i s approximated using Borishansky's c o r r e l a t i o n 5 " for nucleate b o i l i n g forced convection and modified by an e f f i c i e n c y factor to r e f l e c t the 10 geometry. The temperature d i s t r i b u t i o n i s then worked out. in the v e r t i c a l and horizontal sections using carbon conducting paper, i . e . a 2-D technique. For the v e r t i c a l section through the cooling channel, the heat flux estimate i s used as the boundary condition, but thi s i s raised 30% for the horizontal section a r b i t r a r i l y . At the channel the boundary condition i s a fixed temperature condition. The maximum temperatures calculated are then compared to the maximum temperature tolerable and the design accepted on t h i s basis. The design pays less attention to the influence of a x i a l heat flow. The most important parameters are the heat flux and the heat transfer c o e f f i c i e n t to the coolant. In a sample ca l c u l a t i o n the maximum mould face temperature reached 350°C. The tolerable temperature was given as 500°C. No material was mentioned but from the value of the thermal conductivity used i t is assumed that the material was copper. This represents very high mould temperatures when compared with western design. It must be pointed out, however, that most of t h i s work was done more than 10-15 years ago with probably limited computing resources. The approach to evaporative mould design i s sim i l a r . In practice, mould thicknesses are lower (about 5mm) and the material employed i s s t e e l . Such moulds apparently are not used to-day. 11 1.2.4 Water Cooled Moulds Used Elsewhere Water cooled moulds are also employed in other processes. Vacuum arc melting moulds are also designed using the same techniques. However, the heat fluxes are lower and the region of high heat flux i s narrow which, p r a c t i c a l l y , means that lower water flow v e l o c i t i e s are used. There i s creep deformation, but the number of heats obtained i s s u f f i c i e n t to economically amortise the cost of repair and the cost of the mould. 2 1 The same remarks are probably applicable for electron beam melting moulds. Moulds in continuous casting of steel are subject to similar heat fluxes and have been subjected to intensive s t u d y . 2 3 ' 2 4 E a r l i e r , the problem of mould d i s t o r t i o n was attempted to be solved by lowering temperatures in the mould w a l l . " 1 Attempts are now being made to reduce the p o s s i b i l i t y of intermittent b o i l i n g and to gauge the effect of thermal s t r e s s e s 2 4 ; and so design to reduce product r e j e c t i o n . 1.2.5 Probable Mould Type In ESC In e l e c t r o s l a g casting, the casting i s not expected to be very t a l l . Further there w i l l probably be frequent changes in section. One piece moulds, i f employed, w i l l have to be cast. Such moulds may have casting defects ( e.g. blow-holes ) and i t w i l l be d i f f i c u l t to either jacket them or to have the cooling channels i n t e g r a l l y cast. Therefore, i t i s l i k e l y that the mould w i l l be made of short sections joined together by 12 mechanical means. This system favors channel cooling of moulds since i t i s thought that water-jacketing of each section w i l l be too expensive. It i s unlikely that electroslag casting production runs w i l l be very large meaning that copper moulds would probably make this process noncompetitive. The only other material, commonly available, with a s u f f i c i e n t l y high thermal conductivity i s aluminium. This change w i l l cause higher temperatures in the mould w a l l 2 5 , but there i s no report of any untoward e f f e c t on the process or the cast metal. 1.3 Surface Quality Electroslag ingots have a good surface and t h i s was passed off i n t u i t i v e l y as being due to a thin slag skin. M i t c h e l l and E t i e n n e 2 6 suggested that, as in other s o l i d i f i c a t i o n processes, t h i s was due to the presence of a c y l i n d r i c a l l i q u i d metal head. This i s generally, accepted with variations , primarily on the existence or non-existence of a l i q u i d f i l m between the slag skin and the s o l i d i f y i n g metal. A bad surface can come about only as a res u l t of meniscus s o l i d i f i c a t i o n , leading to intermittent r i p p l i n g . 13 1.3.1 Slag Skin C r y s t a l l i s a t i o n The use of water cooled moulds makes the formation of a s o l i d slag layer, c a l l e d the slag skin, on the mould surface inevitable. There has been l i t t l e published information on the formation of slag skins by the s o l i d i f i c a t i o n of d i f f e r e n t slags. Sharapov et a l , 6 0 examined the role of the slag skin. They divided slags into two major groups, those with two components, and those with many components. Both types of slag caused three layers to form in the slag skin with the two layers nearest the metal held to form on reheating of the o r i g i n a l skin. The t h i r d layer closest to the metal was said to be l i q u i d t i l l the metal s o l i d i f i e d . Multi-component slags were preferred in moving mould processes because they permitted a thicker t h i r d layer. Kamensky et a l , 2 7 considered the thermal s t a b i l i t y of the slag and the structure of i t s boundary, the slag skin. The thickness was governed by heat transfer. On petrographic examination, of a 70:30 slag skin they found fiv e layers. The f i r s t layer had calcium aluminate c r y s t a l s with the gaps f i l l e d with f l u o r i t e with the former as high as 70%. The second layer has columnar corundum c r y s t a l s (upto 80%) and the t h i r d layer has f l u o r i t e (92%, rest glass) oriented in the di r e c t i o n of heat flow. Layer 4 i s the same as layer 3 with a small quantity of needles of corundum. Layer 5 has extremely long and large corundum needles (12mm in length) oriented in the heat flow d i r e c t i o n . These needles were much longer than those to be found in granulated slag. They observed that t h i s structure 14 meant that the phenomenon of f r a c t i o n a l c r y s t a l l i s a t i o n was ocurring. The higher alumina content.of the skin ( twice as high ) was attributed to thermodiffusion of fluoride ions and enrichment by ions of the type ( A l „ 0 7 ) - 2 . Korousic and O s t e r c 2 8 also made a mineralogic analysis. They found that the phases c r y s t a l l i s i n g out did not follow the phase diagram. The slag skin was stated to form under impeded metal ion exchange with consequent unbalanced f r a c t i o n a l c r y s t a l l i s a t i o n . It is unfortunate that they f a i l e d to study v e r t i c a l sections close to the metal meniscus. Medovar et a l , 1 mention three layers with no abrupt boundaries. The la s t layer i s said to have a near eutectic composition. The f i r s t layer has an equilibrium composition and forms at very high cooling rates. The middle layer contains non equilibrium proportions of phases and multicomponent slags are more prone to t h i s . B e l l and M i t c h e l l 2 9 also observed the slag skin of slags in the calcium fluoride-alumina system. The slag skin was said to be a near equilibrium phase and they deduced the di s s o l u t i o n and melting of the skin, as the metal meniscus approached, from X-ray maps. M i t c h e l l 3 0 has analysed the slag skins of various slags using X-ray d i f f r a c t i o n . In a 70:30 slag skin CaF 2, A l 2 0 3 and CA6 were found. A 40:30:30 slag caused C 1 2A 7 and C^k7F to c r y s t a l l i s e and a 55:35:10 slag skin had CaF 2 and CA 2. 15 1.3.2 Slag Skin E f f e c t s On Surface Quality In the West, the suggestion of M i t c h e l l and E t i e n n e 2 6 has been generally accepted. The slag composition i s important and i f a high melting phase prec i p i t a t e s out during s o l i d i f i c a t i o n then such a slag i s conducive to the formation of a good surface. If the l i q u i d head of metal i s present, then i t i s held that some of the p r e c i p i t a t e d high melting phase i s redissolved into the slag leaving a smooth surface for the l i q u i d metal to s o l i d i f y against. The Russian workers 1 v i s u a l i s e a thin l i q u i d f i l m between the l i q u i d metal and the s o l i d i f i e d slag skin. They envisage some remelting and/or dissolution of the slag skin and rippled surfaces are attributed to meniscus s o l i d i f i c a t i o n . Some Japanese researchers' o p i n i o n s 3 1 tend to agree with the Russian view, others 9 favor the M i t c h e l l viewpoint. The surface of the ingot i s said to be controlled by the freezing point of the a l l o y and the nature of the slag. A high melting phase p r e c i p i t a t i n g out favors the formation of a smooth skin. Speculation in industry suggests that smooth skins may be encouraged by running the mould warmer. However, M e l l b e r g 3 2 has pointed out that the heat flux in the molten head region of the ingot i s very high and t h i s may act to control the superheat in the l i q u i d . This would mean that running a mould warmer may have no e f f e c t . In i n d u s t r i a l practice, the problem of a poor surface has been generally dealt with by increasing the power input, but th i s i s limited due to adverse effects on s o l i d i f i c a t i o n 16 structure. It i s also known that a higher f i l l - r a t i o favors a better surface in fixed mould practi c e . An increase of slag depth i s also favorable but t h i s , i t i s stated, i s unexpected. 3 3 The mould, i t appears, has no e f f e c t . Shevtsov et a l , 3 4 suggest that a mould can have no eff e c t on slag skin formation because no mould can be made that has a face temperature in excess of the slag melting point based on the premise that the slag skin thickness adjusts i t s e l f to match the heat flux from the bulk of the bath to the mould. 1.4 Objectives The foregoing discussion suggests strongly that moulds to be used in electroslag casting w i l l consist of channel-cooled sections made from aluminium. It i s noted that electrodes to be used w i l l normally be round bars. With t h i s in mind, investigations were ca r r i e d out to determine i f such moulds could be designed to ensure smooth surfaces in electroslag casting. It was of interest to design moulds to have a given l i m i t i n g temperature at the mould surface. Besides the obvious necessity of keeping the mould cool enough to prevent disastrous melting, there were suggestions that temperature i s the c r i t e r i o n for mould design. It was intended to determine i f th i s c r i t e r i o n was correct or i f there were other more important c r i t e r i a in mould design. 17 I I . EXPERIMENTAL This chapter i s divided into two sections. The f i r s t deals with experiments made on the 1 tonne electroslag furnace and the second i s concerned with observations of the slag skin and the slag cap as obtained after melting. 2.1 Experiments On The Furnace It was important to measure temperatures in the mould wall at d i f f e r e n t locations, since data was needed to compare with model predictions. A heat flux sensor also had to be developed since i t was known that the heat flux varied with position in non-symmetrical geometries and values of heat flux were needed as input to mathematical models of heat transfer in the mould. 2.1.1 The Electroslag Furnace At UBC A d e t a i l e d description of the furnace has been presented e a r l i e r . 3 5 ' 2 5 The furnace i s capable of melting ingots and castings up to 1 tonne. There are two transformers, the f i r s t rated 250 KVA with voltage dropping from 12.5KV to 600V and the second d e l i v e r i n g e l e c t r i c a l power at between 25V to 60V in steps of 2.5V. The electrode holder i s water cooled and moves v e r t i c a l l y up and down between aluminium r a i l guides. The holder carriage i s suspended from a chain, the chain driven by a motor whose 18 speed i s continuously variable up to 163mm/min. The melting rate was controlled by varying the electrode descent v e l o c i t y . The. current was i n d i r e c t l y controlled by the electrode v e l o c i t y and a saturable core reactor was not used. V e r t i c a l channels were d r i l l e d in the mould and horizontal holes connected the channels to brass nipples at the ends of the mould. The open holes at the extremities of the channel were sealed by welding. The moulds were mechanically assembled and cooling water brought to the moulds through rubber hoses with thermometers monitoring the i n l e t and outlet temperature. Both aluminium and copper sections were used in the course of the work. A picture of the assembly i s shown in f i g . 3 and 4. Sketches of the mould sections in which the thermocouples were embedded are given in f i g s . 5(a,b,c) and 6. 2.1.2 Temperature Measurements Thermocouples are the most convenient means of temperature measurement in moulds and have been used extensively by previous workers. The maximum temperature expected to be measured was about 300°C. Hence, copper-constantan thermocouples were chosen. These were available in thicknesses of 10 thou and sheathed in f i b r e g l a s s , which also served as insulation with the diameter of the assembly being about 1mm. Several holes of diameter 1.25 mm were d r i l l e d in a mould section as shown in figures 5(b,c) and 6. The type A thermocouples were expected to give an idea of the temperature 19 Figure 3 - Sketch of Mould Assembly Figure 4 - Views of Mould Assembly 20 a Hot Foce Type A X Y Pos. No. (in cms.) 1 3.57 95 3 4.61 .95 5 3.11 1.43 7 3.58 1.91 9 460 1.51 11 4.66 2.05 Typ« A Z X Pos. No. cm cm 1 11-46 3 57 3 12 67 461 5 13 97 311 7 1522 356 9 1652 4.10 11 17 78 4-66 o Type B Figure b 5 -c Dimensions of Mould Section used 21 13-5 1 20 1 20 10 dimensions in mms. section height = 254 mm T/C No. X ( y in mm z ) U 5-1 9-5 101-6 15 5-8 3 6 127 16 3-5 3-6 1524 COPPER SECTION USED Figure 6 - Dimensions of Copper Section d i s t r i b u t i o n in the mould, while type B thermocouples were to confirm the assumption of steady state heat flow, i . e . type B thermocouples were expected to give the same temperature h i s t o r y under steady state conditions. The thermocouples were welded outside and soldered to the bottom of the holes using STRONGSET No. 509 solder. The cold junctions were encapsulated in insula t i o n tape and placed in test-tubes which had some mercury in them to promote heat transfer and were kept at 0°C by p a r t i a l immersion in an i c e -water mixture. The measuring junctions were connected to a Texas Instruments FM6WB multi-channel potentiometer-recorder. 22 2.1.3 The Heat Flow Sensor Considerable d i f f i c u l t y was found in f i t t i n g published values of heat flux to the measured temperatures. It was evident that heat flow was d e f i n i t e l y three dimensional and t h i s created the need to develop a heat flow sensor. The heat flow sensor had to be thermally insulated from the rest of the mould and had to be as thin as p r a c t i c a l to permit accurate heat flux measurement. Temperature measurement of water at the i n l e t and outlet was considered and rejected since i t required too great a s e n s i t i v i t y in available measuring instruments. The heat flux sensor was f i n a l l y designed as shown in figure 7. The thermocouples were copper-constantan and the sensor was heavily water cooled so that no more than 0.5 degrees Celsius change in temperature was expected between the i n l e t water and the water at the outlet. The sensor was wrapped t i g h t l y in FIBREFRAX to thermally insulate i t . Such sensors have been widely used by Russian researchers. 1 It i s noted that the metal probably constitutes a neg l i g i b l e thermal resistance and so the choice of material for the sensor may not matter. 2.1.4 V e r i f i c a t i o n Of Sensor Data The temperature difference between the thermocouples was expected to be d i r e c t l y proportional to the heat flux. Published values of the heat transfer c o e f f i c i e n t suggested that the Biot number was about .02 to .04. Here i t i s pointed out that the Biot number i s normally used in cases where there i s no 23 o G - _ 50 mm-— — » 0 . . . Thermocouples Heat Flux Sensor Figure 7 - Sketch of Heat Flux Sensor 100 80 60 I 20' 0 0 10 20 30 40 50 60 70 Diff. in Temperature. *C Figure 8 - For V e r i f i c a t i o n of Sensor 24 sink as i s present in the sensor, but an extension of the idea may give an insight into the effectiveness of the sensor, approximate as i t may be. It i s suggested that for a straight l i n e v a r i a t i o n of temperature to exist in the sensor the Biot number at the source should be small and that at the sink should be large. Since the heat transfer c o e f f i c i e n t at the coolant interface i s about 4 orders of magnitude higher than that at the hot face, i t appears l i k e l y that a straight l i n e v a r i a t i o n of temperature w i l l occur. The response time, L 2/A, where L i s the c h a r a c t e r i s t i c length and A the thermal d i f f u s i v i t y , was calculated and found to be about 3 s. The rate of r i s e of ingot i s about .05 cm/s which suggests that the sensor w i l l probably respond quickly enough. To v e r i f y t h i s estimate the difference between the two temperatures (proportional to the heat flux) was plotted against the temperatures. Provided that the sensor was good the plot should be a straight l i n e s t a r t i n g at the temperature of the coolant at the i n l e t . This i s shown in figure 8. It was considered that the data was adequate to be used after i t had been normalised using figure 8. 2.1.5 Experimental Results A summary of the 15 melts made i s presented in table I I . Temperatures were measured at d i f f e r e n t positions in the mould. Later, when i t became apparent that the temperature could be predicted reasonably well, only the temperature of position 3 and 11 were recorded. Run 8 had the mould section reversed so TABLE I I Run # CaF„ wt% S CaO wt% l a g 8 kg. A l 0 wt % SIO wt % MgO wt % Electrode Dia. (mm) S p e c i a l Features 1 70 — 30 _ _ 75 T r i a l 2 70 — 30 _ _ 75 Mould Temperatures and Heat Flu x Measured 3 70 15 15 _ _ 113 Thicker e l e c t r o d e , temp., heat f l u x measured 4 70 _ 30 _ •_ 113 I I 5 50 20 30 113 M 6 49 16 17 12 6 113 Slag leak, melt aborted 7 49 16 17 12 6 113 Thick e l e c t r o d e : temp., heat f l u x measured 8 70 — 30 _ _ 113 Mould s e c t i o n reversed, power f a c t o r = 0.96 9 70 — 30 _ _ 113 Melt aborted. Cable overheated 10 70 — 30 — _ 113 Copper s e c t i o n replaces aluminum 11 70 - 30 - - 113 Run #10 repeated to avoid i n f l u e n c e of power f l u c t u a t i o n s 12 70 — 30 — — 113 Grooved s e c t i o n used w i t h T/Cs and f l u x sensor 13 70 — 30 - — 113 F l u x sensor at corner of box 14 70 — 30 — — 113 F l u x sensor at corner of box 15 70 - 30 - - 113 E l e c t r o d e as i n F i g . 9 26 - Holding Stub 50 mm f 400 mm Starter Stub Figure 9 - Electrode used in Run 10 RUN 3(exp.| Figure 10 - Temp. History of RUN 3 27 £ s. 0 8 16 24 32 40 48 56 Time. mins. Figure 11 - Smoothened Temp. History of Run 3 94=—, 1 , , . . . , , . , , , r-0 8 16 24 32 40 48 56 Time. mins. Figure 12 - Calculated Temp. for Run 3 28 Figure 13 - Temp.History of Run 2 RUN 2lcalc> 0 8 16 21 32 40 48 56 Time, mins Figure 14 - Calculated Temp. for Run 2 29 Figure 15 - Temp.History of Run 6 Figure 16 - Calculated Temp. for Run 6 RUN lOlexpl 14 • 15 — 16 • 0. !0 20 30 40 50 60 Time. mins. Figure 17 - Temp.History of Run 10 gure 18 - Calculated Temp. for Run 10 31 Figure 19 - Temp.History of Run 8 Figure 20 - Heat Flux Measurement-Run 2 32 0 I . . , r-0 1 2 3 4 Heat Flux. MW.m"2 Figure 21 - Heat Flux Measurement-Run Figure 22 - Heat Flux Measurement-Run 33 Time (mins) i i r Heot F l u » . MW.m" Figure 23 - Heat Flux Measurement-Run 6 RUN II 0 1 2 3 4 Heot Flux. MW.m"2 Figure 24 - Heat Flux Measurement-Run 11 34 0 8 16 24 32 40 48 56 Time. mins. Figure 25 - Sample comparison between Calculated temperature and measured that thermocouples previously placed higher in the mould were now at a lower p o s i t i o n . In this run the power factor was determined to be 0.96. Run 12 was made with a grooved face. Runs 13 and 14 were made to measure the heat flux at the corner of the box section. Run 15 had an electrode as shown in f i g . 9. Figures 10,11,13,15,17,19 d e t a i l the temperature history obtained. The curves have been moved on the plots to simulate the temperature f i e l d s on horizontal planes in the mould. Fig.10 shows the variation obtained with heavy power flu c t u a t i o n . In the other runs the flu c t u a t i o n was not so high, but the curves have been smoothened to reduce the influence of power fluctuation and f a c i l i t a t e better comparision with 35 computed r e s u l t s . It i s clear from runs 2 and 4 that thinner electrodes have lower melting rates and cause lower heat fluxes to the mould than thicker electrodes. The temperature in the mould rose and f e l l smoothly as the slag passed across the measuring plane, except when there was a power i n s t a b i l i t y . The heat flux curve showed two peaks, one at the sl a g - a i r interface and the other at the l i q u i d metal head l e v e l with the peak flux l e v e l at the l i q u i d metal head being much higher than previously reported. (see f i g s . 20 to 24 ) The heat flux at the corner of the square section was determined with some d i f f i c u l t y . In the f i r s t attempt, the flux was found to be much higher than expected from mould temperature measurements. In the second attempt, special e f f o r t s were directed at ensuring that there was no slag penetration around the sensor. It was thought that such penetration would prevent contraction of the slag skin away from the sensor. Apparently, the a i r gap between the slag and the mould played a major role in c o n t r o l l i n g the heat fl u x . Therefore the values of heat flux determined at t h i s position may not be very accurate. What i s clear i s that the average heat flux and the magnitude of the peak flux are both lower. Based on these results the t o t a l heat flux to the mould was calculated and compared to the power input and the match (shown in table III) was found to be sa t i s f a c t o r y . Run 8 was made to v e r i f y that there was some transient heat transfer. The maximum temperature reached by position 3 was found to be higher than that reached by position 11, whereas i t was in the opposite sense in e a r l i e r runs thus demonstrating the 36 Table III - Comparison: Heat flow and Power input Run# Area Rate of ingot r i s e cm/min. Heat flow (int.) kW Power input kW Mould ht. considered cm cm 2 3 4 5 6 248.7 365 215 194 236.8 .37825 .22325 .38372 .40959 .28069 127.8 1 10.7 112.1 107.98 90.02 135.0 1 35.0 1 35.0 135.0 126.0 17.02 13.17 16.11 13.92 1 5.43 existence of transient heat transfer. 2.2 Observations On Slag Skins The formation of the slag skin plays an important part in co n t r o l l i n g heat transfer and in the formation of the metal surface. It was necessary to examine the process of skin formation in order to ascertain i f i t could be modified so that a smooth casting surface could be ensured. With t h i s view, the slag skins formed with various slags were examined. The slag cap of 70:30 slag was also examined. This slag was chosen since i t s phase diagram i s r e l a t i v e l y well known. 2.2.1 Microstructure Some micrographs of the slag skins are shown in figures 26 to 35 from which i t i s clear that the structure consists of three layers. The f i r s t seems to be a mixture being dark in some areas and l i g h t in others. There does not appear to be any pa r t i c u l a r regular structure or any predominant component in the 37 Figure 27 - Spike on 70:30 skin 38 Figure 29 - 50:30:20 slag skin Figure 31 - Skins of complex slag 40 Figure 32 - Industrial 70:30 slag skin Figure 33 - Industrial 50:30:20 slag skin 41 Figure 35 - 70:30 slag cap near top 42 layer. The second, or the middle, layer displays r e l a t i v e l y thick c r y s t a l s oriented toward the ingot casting. There are several voids but most of the space in between the thick c r y s t a l s i s f i l l e d with a darker phase. The c r y s t a l s appear to be the primary phase while the darker phase i s probably the eutectic l i q u i d which s o l i d i f i e d l a t e r . Though the word ' c r y s t a l ' has been used i t i s pertinent to note that in some of the complex slags the middle layer actually looks l i k e a glassy porous mass. The t h i r d layer, the one closest to the ingot casting does have some primary c r y s t a l s , but s u r p r i s i n g l y , has the physical appearance of the secondary phase. This layer was d i f f i c u l t to retain in the sample during polishing and caused considerable d i f f i c u l t y during that procedure. Apparently i t i s only l i g h t l y held against the skin. In the p a r t i c u l a r case of 70:30 slag, there were some spikes of glassy transparent material s t i c k i n g to the slag skin at random i n t e r v a l s . In the other slags especially the 70:15:15 slag, the t h i r d layer was very thin and at some places not v i s i b l e . 2.2.2 Observations On The SEM Following the work under the l i g h t microscope, the EDX detector on the scanning electron microscope was used to i d e n t i f y the elements present in the three d i f f e r e n t layers. Quite obviously d i f f e r e n t slags would present d i f f e r e n t results, but the general pattern became clearer as seen in f i g s . 36 to 43 39. In the case of 70:30 slag the outer l'ayer was inconsistent; in some areas Al would be present and in other areas absent. This was not the case with Ca; i t was always present and the Ca peak was higher than the Al peak. In the middle layer, there was a strong peak for A l , even larger than that for Ca. In the inner layer, close to the ingot, Al was noticeably absent except in the c r y s t a l s that pierced their way through from the middle layer. The spikes of glassy material showed an even more surprising composition; apart from the presence of Al there was Si in s i g n i f i c a n t content and noticeable quantities of S,Fe and Mn. This was surprising because neither the remelting theory nor the l i q u i d f i l m theory could explain the composition of the t h i r d layer or that of the spikes. The structure of the 70:15:15 slag was d i f f e r e n t . The outer layer showed A l , Ca and S i . The middle layer had no Si and the peak due to Al had increased in height. The inner layer, closest to the ingot, showed an inconsistent presence of Al and Si e s p e c i a l l y at the very edge of the skin. There were no spikes on t h i s slag skin. The 50:30:20 slag skin had in i t s outer layer Al and Ca with the Ca peak stronger. The middle layer had two distinguishable phases, one darker than the other. Both Al and Ca were present with the Ca peak s l i g h t l y larger in the dark phase but smaller in the l i g h t phase. The inner layer showed very l i t t l e Al but S i , S and Mn were in evidence. An i n d u s t r i a l slag of 70:30 type used to melt stain l e s s 44 ••••'X • */...• • • - • -'".v. -f outer layer ' . "•' **' middle layer ' ' V.Vt/'.V'v'.' • ^•?v>V.'^ v-cvv.. ... o .': . . . inner layer o> • • . : : A. • • * • ' V .** . • • • base of spike • . .'. • • • .v. •. •/ spike .-. .. .••»."'•"* A| Si Ca Ca Mn Fe '' t f t f 'k Figure 36 - EDX analysis of 70:30 skin 45 •A o u o outer layer • A . •• • middle . 4 2 J> • •. inner layer Al Si Ca Ca II LL Figure 37 - EDX analysis of 70:15:15 skin 46 \~ \~y % o o • tn _ i - •• •* inner^ layer' light phase (middle) •V -.V*.- .'• dark phase (middle) outer layer Mg Al Si S Ca Ca Mn I L L ! tt T Figure 38 - EDX analysis of 50:30:20 skin 47 . . " outer layer ' ' ' * -V^'rvVV/V *\ I light phase (middle) dark phase (middle) .'• -.\Vv. V:v^. V - .A. - . . . . . inner layer * rv.V.,^;..^**. Mq Al Si Ca Ca Mn Fe f t t t t t t Figure 39 - EDX analysis of complex slag skin 48 steel was examined and displayed much the same pattern with Cr and T i appearing on the inner layers r e f l e c t i n g the effect of the metal melted. A slag having 12% S i 0 2 , 6% MgO, 16% A l 2 0 3 , 17% Cao, and the rest CaF 2 showed the presence of Mg, S i , A l , Ca and Mn, (the last a pickup from the metal) in the outer layer. The middle layer had two phases, the darker of which had A l , S i , and Ca ; the l i g h t e r had Mg, Al traces of Ca and small amounts of Mn and Fe. The inner layer had a l l the elements mentioned though the peak of Ca was dominant. 2.2.3 The 70:30 Slag Cap Near The Liquid Metal Meniscus It was apparent that the middle layer consisted mainly of primary c r y s t a l s . It was also clear that these c r y s t a l s were Table IV - Chemical analysis of some skins Sample CaF 2 Cao A1 20 3 wt% wt% wt% 70:30(bulk) 80.7 — 15 70:30 skin 60 to 70 — 30 to 40 70:15:15 skin 38.4 17.9 17.5 50:30:20 skin 43.4 29.7 29.9 much thicker than the c r y s t a l s that would have normally precipitated (those observed in the slag bulk) and d e f i n i t e macrosegregation was observed. Micrographs were taken and the proportion of primary phase in the skin and in the slag bulk determined. The proportion of primary phase was 4-5 times 49 higher in the skin than in the bulk. Further, there was a clean break at the skin and no evidence that remelting or dissolution took place at the meniscus. The middle layer started growing about 1cm below the slag top and increased in thickness t i l l the metal meniscus in a nearly linear fashion. An e f f o r t was made to quantit a t i v e l y determine the proportion of elements present in the various phases using electron probe microanalysis. The presence of oxygen and fluorine meant that the use of a l i g h t element spectrometer was necessary. The f r a g i l e nature of the specimen coupled with the presence of voids in the specimen probably reduced the s e n s i t i v i t y so that r e l i a b l e analysis became d i f f i c u l t . Analysis by atomic absorption was possible and the analysis i s given in Table IV. Apparently much of the alumina did not dissolve in the slag and the composition of the 70:30 slag i s actu a l l y close to 80:20, but the increase of alumina in the skin i s evident. The l a t t e r analysis cannot, of course, indicate the individual composition of the three layers. 50 I I I . CALCULATIONS Several simple models were formulated in an e f f o r t to explain the observations. The f i n i t e difference technique was used, with some success, to calculate mould temperatures. In a speculative c a l c u l a t i o n to determine the effect of grooves in the mould face on heat transfer to the mould, the f i n i t e element method was used. 3.1 Calculations In 2 Dimensions The f i n i t e difference technique i s now a standard method to solve the d i f f e r e n t i a l equation for steady state heat transfer. Good explanations are given in many textbooks, e.g. those by Kreith and B l a c k , 3 8 by Patankar 3 6 and by Carnahan, Luther and W i l k e s . 3 7 I n i t i a l l y , a rough estimate of the temperature f i e l d was needed to determine the positioning of thermocouples in the mould. Therefore, a 2-D steady state heat flow s i t u a t i o n in a. horizontal s l i c e of a f l a t mould section was examined. The f i n i t e difference d i s c r e t i s a t i o n was made, solved and the results are shown in figure 40 (for a heat flux of 25 cals/cm 2K). It i s clear that the temperature gradients near the hot surface are very high and any error in placement of a thermocouple could lead to large errors in measurement, making them unreliable. It i s also seen that the region marked B l i e s in a area of low temperature gradient which means that i f one d r i l l s a hole from the rear face as shown there i s l i t t l e chance 51 ure 40 - Temp. f i e l d c a l c . with 2-D assumptions Hot Face Using symmetry, only portion bounded by heavy lines discretised. Arrows show heat flow. Figure 41 - The s l i c e which was d i s c r e t i s e d 52 of inaccuracy in temperature measurement in the region. Type B positions, i t seems, may be used for v e r i f i c a t i o n . Measured values of heat flux were used in the c a l c u l a t i o n and i t was quickly discovered that a 2-D model was inadequate since i t predicted extremely high temperatures at the hot face. It was evident that computation had to be extended to three dimensions. 3.2 Calculations In 3 Dimensions The next step was to determine i f a quasi steady state model could calculate temperatures in the mould that were reasonably representative. The results of the previous section suggests that an accuracy of no more than about 5 Celsius degrees in temperature measurement may be expected. Following the type of assessment used for the heat flux sensor a similar l i n e of reasoning indicated that a mould section such as considered was not very suitable for the assumption of quasi steady state heat transfer. Calculations of a "Biot" number indicated that the condition of quasi steady state may exist only in a body of size no greater than 4-5 cms. The thickness of the mould was 5cms and the height about 25 cms and so there i s probably an inherent error in making quasi steady state assumptions. However, there i s a r e l a t i v e l y slow r i s e of the slag across the mould face. Further other workers have shown that the heat flux to the mould in the slag region i s r e l a t i v e l y constant 1 and i t i s possible that the maximum 53 temperatures computed on the mould face may not be inaccurate, i . e . a steady state c a l c u l a t i o n may predict maximum temperatures with acceptable accuracy. It i s also to be noted that three dimensional unsteady state computing i s p o t e n t i a l l y extremely costly and hence unsuitable in the f i e l d . A three dimensional quasi steady state c a l c u l a t i o n was made. The technique used matrix solutions as well as i t e r a t i o n in an e f f o r t to balance computing costs against memory requirements. A single run took about 20 seconds CPU time on the AMDAHL computer at the University of B r i t i s h Columbia. A program l i s t i n g in FORTRAN i s given in Appendix C. The results of the program written for 2-D approximation were used as the i n i t i a l solution estimate for the program written assuming 3-D steady state heat flow. The input required i s the heat transfer c o e f f i c i e n t to the coolant, the heat flux d i s t r i b u t i o n , the thermal conductivity of the mould, the physical dimensions, the coolant flow rate, and the coolant temperature. In the ca l c u l a t i o n , the heat transfer c o e f f i c i e n t to the cooling water was estimated by both the Sieder-Tate 5 5 c o r r e l a t i o n as well as the r e l a t i o n given by Petukhov. 5 6(App.B) There was a small difference but i t was seen l a t e r that the mould face temperatures were not s i g n i f i c a n t l y changed by these va r i a t i o n s . The heat flux into the mould face was taken from the experimental data. Iteration to the nearest degree was deemed s u f f i c i e n t since the accuracy of measurement was no greater than that. The thermal conductivity was assumed constant and equal to the value at 200°C for aluminium but for 54 copper a lin e a r v a r i a t i o n with temperature was permitted. Since r e l a t i v e l y good agreement with experimental data was observed, the program was not refined to allow a better estimate of the thermal conductivity. The re s u l t s are given in figures 12, 14, 16, 18 (p.27-30). The temperatures were calculated using 3 s p a t i a l dimensions but the v e r t i c a l dimension (z) can be converted to time since the rate of r i s e of the ingot i s r e l a t i v e l y constant. The results have been given with respect to time to f a c i l i t a t e comparison. It i s to be noted that points within the mould in the same position with reference to the slag have approximately the same temperature hist o r y . This i s true except for the 2-3 cms at the base and at the top of the mould. The temperatures calculated are not quite the same as the measured temperatures (as expected) but are close enough to be considered reasonable and th i s comparison i s v i s u a l l y available in figures 11 to 18 and 23. The temperature gradient with time is also reasonably correct. It i s f e l t that such calculations may be used to predict maximum temperatures in mould sections though there i s l i t t l e doubt that transient heat transfer has a discernable effect on temperature f i e l d s . This i s confirmed by the fact that steady state calculations under-predict temperatures at type B positions and also by the res u l t s of run 8. Following the acceptable success of the model the effect of variation of mould material, channel positioning and channel spacing were made. As expected, the use of copper lowered hot face temperatures dramatically. This i s seen in f i g . 17 and 18 55 (p.30). Moving the channels closer reduced the face temperature and made the temperature gradient more uniform (fig.43 a,b). If the channel i s brought closer to the hot face, the face temperature is reduced but the gradient increases and becomes less uniform (fig.43 c,d). There i s also a high temperature at the cooling channel and t h i s may create b o i l i n g thus changing the cooling conditions. A small change (10%) of the heat transfer c o e f f i c i e n t to the coolant brings only small changes in temperature d i s t r i b u t i o n suggesting that the approximation used to estimate the c o e f f i c i e n t i s acceptable (fig.42 b,c). It i s to be noted that the heat transfer c o e f f i c i e n t depends on the mean temperature of the coolant f i l m according to the method used and so should vary over the length of the channel. A c a l c u l a t i o n with an a r b i t r a r y variation was also performed and showed i n s i g n i f i c a n t changes (fig.42 d). Variation of i n l e t water temperature also had only a small e f f e c t , but i t i s noted that i f the temperature r i s e s above 25-30°C at the same cooling v e l o c i t i e s , then there i s a strong p o s s i b i l i t y of b o i l i n g occurring (fig.44 a,b). The horizontal plane chosen is that where the hot face reaches i t s highest temperature. c) h raised 10% d) a r b i t r a r y h Figure 42 - Calculated Temp. f i e l d s 57 a,b) reduced inter-channel distance c,d) channel closer to hot face Figure 43 - Calculated Temp. f i e l d s 58 a) water at 25°C b) water at 10°C Figure 44 - Calculated Temp. f i e l d s sketch of section token for finite element calculations Figure 45 - Section taken for F.E. c a l c u l a t i o n 59 90° V-groove 304 (cals/cm^s) 10-+ — — , , r -0 10 20 30 40 °/o Relative Surface Heat Flux tcals/cm^s) 0 U - groove 10 20 30 %> Relative Surface LO Figure 46 - Calculated results of speculation 60 3.3 Speculation On The Eff e c t Of Grooves On Surface Quality The assumption of steady state heat transfer coupled with knowledge of the slag skin thickness and i t s thermal conductivity indicated that the a i r gap had a smaller thermal resistance, though not i n s i g n i f i c a n t , compared to that of the slag skin (see Table I, VI and VII the l a t t e r on p. 67-68). Thus, a s l i g h t increase of the a i r gap could reduce the heat flux s i g n i f i c a n t l y . It i s known that rough mould surfaces in metal ingot casting lead to smoother surfaces, presumably due to lower heat fluxes. It was thought that machining grooves in the mould surface might have the same e f f e c t . There were two assumptions made. The f i r s t was that the slag would not enter the grooves since i t had been observed that slag did not enter the space between two mould sections, which could be as large as 2mms. The other was that the slag skin thickness would not change as the molten metal head approached. To determine the eff e c t of groove spacing on heat flux, a f i n i t e element model of a s l i c e of mould was formulated, a sketch of which i s shown in figure 45. The FORTRAN l i s t i n g of the program i s given in Appendix D. The results of the speculation are given in figure 46. It seemed that the heat flux could be e f f e c t i v e l y reduced by up to 25% i f the grooves covered 50% of the hot face. Apparently, there was need to make experimental investigations since i t seemed possible to design a mould to permit a molten head of metal, to exist at lower power input than hitherto and t h i s was the basis for run 12. The experiment was unsuccessful since i t was discovered 61 that to prevent s l a g from e n t e r i n g the grooves , i t i s necessary to quench the s l a g extremely r a p i d l y . The genera l h e a t i n g up of the mould w i t h the r i s e of the s l a g makes t h i s i m p o s s i b l e . No exper iments w i t h a groove width s m a l l e r than 1 mm were attempted s i n c e i t was thought that such f i n e grooves would be expensive to machine and so would be i m p r a c t i c a l . 62 IV. DISCUSSION 4.1 Failu r e Of 2-D Model Russian workers 1 have calculated temperature f i e l d s using carbon conducting paper. This presupposes that heat flow can be adequately described as two dimensional. While t h i s may be permissible for a v e r t i c a l section exactly between two channels i t i s not s u f f i c i e n t for a horizontal section. This was re a l i s e d and i t was suggested by these workers that the heat flux be a r b i t r a r i l y raised by 30%. No other predictions of temperature f i e l d s are available for channel cooled moulds. Channel cooled moulds are used in continuous casting. For th i s case there are some l i t e r a t u r e studies. Chizhikov et a l 4 0 and Simonov et a l 3 9 accepted 2-D steady state approximations. Young, Sunaryo and C r o s s 4 1 used temperature measurements in a v e r t i c a l section exactly between two channels as boundary conditions for the modelling of horizontal sections. This was probably done because of li m i t e d computing resources. Alberny et a l 4 2 assumed n e g l i g i b l e v e r t i c a l heat flow in their measurements of heat flux. V e r t i c a l heat flow has often been neglected in the past in measurements of heat flux in both the electroslag process and in continuous casting. This approach i s flawed since appreciable heat flux gradients e x i s t . The extent of error can be examined by considering a rectangle in the X-Y plane with boundary 63 conditions as shown in f i g . 47. The error i s negligible in the central region of the rectangle but varies from 0 to 10% in the extreme regions. A hypothetical heat flux v a r i a t i o n was assumed (fig.48) and cal c u l a t i o n s performed to check the error. The results are shown in Table V. It is clear that i f the thermocouples used to measure heat flux are placed close to the hot face the error i s low (about 10%) and within experimental error. However, i t i s l i k e l y that the sharp peaks observed in th i s study would not have been uncovered since the type of error associated with the measurements would have obliterated them. It i s also clear the heat flux in the X-direction reduces with increasing thickness f a l l i n g about 15-20% from the hot face value. This indicates why temperature f i e l d s must be calculated in 2 dimensions even i f heat flux measurements are made using 1-D approximation. In the case of channel cooled moulds, the reduction of symmetry necessitates an extension to three dimensions. In the early stages of t h i s work, i t was seen that 2-D approximations of temperature f i e l d s were grossly inaccurate. The heat flux from the slag was an order of magnitude higher than that from the ingot or the region above the slag which caused appreciable v e r t i c a l heat flow and meant that a 2-D approximation would not be s a t i s f a c t o r y . 2-D models of v e r t i c a l sections are acceptable for thin jacketed moulds or i f the cooling channels are extremely clo s e l y spaced, i . e . where the isotherms are f l a t and p a r a l l e l to the hot face. Such moulds w i l l probably not be acceptable in ele c t r o s l a g casting from 64 Figure 47 - Error in heat flux measurement Figure 48 - Hypothetical heat flux v a r i a t i o n 65 Table V - Possible Error in heat flow measurement a. Measuring T/C's 3 and 11 mm from hot face Position Actual Flux 'measured' heat flux in mould of (cm) Cals/cm 2 0C thickness, cm. 1.5 2.0 2.5 4.0 0 5 5.46 5.62 5.74 6.02 6 10 11.5 11.51 1 1 .54 11.61 12 35 32.01 31 .99 31 .93 31.81 17 35 33.63 33.55 33.46 33.26 24 21 21 .08 21.11 21.13 21.18 30 9 1 1 .47 11.51 1 1 .58 1 1 .74 b. 'Measurement' at cold face. 0 5 5.92 6.58 7.33 9.73 6 10 12.86 13.56 14.16 15.16 12 35 29.84 28.87 28.01 25.94 1 7 35 32.59 31.9 31 .21 29.2 24 21 21.16 21 .28 21.41 21 .82 30 9 13.27 14.17 1 4.99 17.12 reasoning given in section 1.2.5 and thus the need to examine 3-D models arose. 4.2 Apparent Success Of 3-D Quasi Steady State. Model Three dimensional approximation i s p a r t i a l l y sucessful, as explained hereunder. The correlations developed for estimating the heat transfer c o e f f i c i e n t to the cooling channel are r e l a t i v e l y accurate only for f u l l y developed flow in the channel and they are r e s t r i c t e d to conditions where the heat flux does not vary substantially -along the channel. When the heat flux does vary the methods developed by S l e i c h e r 4 3 and o t h e r s 4 4 ' 4 5 ' 6 2 for the ca l c u l a t i o n of l o c a l heat transfer c o e f f i c i e n t s should be employed. In 66 theory, even t h i s i s not good enough since the heat flux varies with position even in a horizontal section and in any case the flow i s not f u l l y developed. It i s fortunate that substantial a x i a l conduction occurs, so that the temperatures in the mould near the hot face are not very sensitive to changes in the heat transfer c o e f f i c i e n t . This means that correlations such the Petukhov r e l a t i o n 5 6 may be used. However, i t i s pointed out that b o i l i n g can cause drastic changes and in t h i s case more complex calculations are necessary. The other question regarding the u t i l i t y of a 3-D c a l c u l a t i o n comes from the fact that the heat flow i s a c t u a l l y time dependent. However, the experimental results show that the arguments brought forth e a r l i e r are j u s t i f i e d , to some extent. In i n d u s t r i a l conditions, i t would seem unnecessary to go to the expense of making more complex calculations e s p e c i a l l y i f the motive i s to l i m i t the maximum hot face temperature. The impetus for 3-D unsteady state calculations may come from the the need to calculate thermal stresses. A knowledge of thermal stresses i s necessary since the major cause for mould retirement i s deformation by creep caused by such stresses. However, i t i s noted that the temperatures and temperature gradients may be predicted f a i r l y well using quasi steady state c a l c u l a t i o n s and so i t may be worthwhile to make measurements and compare the results obtained. If the result is precise enough, no refinement of the assumptions may be necessary. The most important variable, over which the designer has l i t t l e c o n trol, i s the heat flux to the mould. There have been 67 many measurements of the heat flux under various conditions, but as yet i t does not seem possible to predict the heat flux for an arbi t r a r y shape. There i s need to develop a method of making r e l i a b l e estimates using simple techniques before such an approach may be used in day to day work. It is suggested that the sharp flux peak observed may be neglected since the peak i s very narrow and the heat flux approximated as a broad band of high heat flux. 4.3 Factors Affecting The Heat Flux To The Mould Table VI - Temp, drop through slag skin a. Heat Flux = 65 cals/cm 2s slag skin Drop in Temperature °C thickness(cm) .006* .008* .010* 0.15 1 625 1218 975 0. 12 1 300 975 780 0.10 1 083 812 650 0.08 866 650 520 b. Heat Flux = 35 cals/cm 2s 0.5 2916 2187 1750 0.4 2333 1750 1400 0.3 1750 1312 1050 0.2 1 1 66 875 700 0.1 583 437 350 * Thermal conductivity in cals./cm.s.°C The thermal resistances to heat flow to the mould from the slag or the ingot are the slag skin and the a i r gap, presumed to exist between the slag skin and the mould. Simple f i t t i n g of 68 Table VII - Heat flow through a i r gap Air Gap (cm) Drop in Temperature, °C for heat flux(Cals/cm 2s) of 65 35 0.1 0.05 65 000 32 500 35 000 17 500 0.001 0.0005 0.0001 650 325 65 350 175 35 A l l units in the CGS system. Thermal Conductivity of a i r taken to be 0.0001 cals./cm.s.°C the known physical properties of a i r and the s o l i d slag to the observed heat fluxes (as in Tables VI and VII) indicates that both the slag skin and the a i r gap are s i g n i f i c a n t resistances. It appears that the a i r gap thickness i s of the order of 10 microns and the temperature drop across i t roughly 400°C the rest of the drop of roughly 700~800oC across the slag skin. It i s pointed out that t h i s i s a very approximate analysis and that the " a i r gap" written about i s only an average value and in fact probably varies dramatically from point to point. In any case one cannot neglect either resistance. The examination of a 70:30 slag cap indicated that the slag skin i s not a simple frozen layer and no remelting or diss o l u t i o n as has been previously stated to occur was observed. The primary c r y s t a l s ( A l 2 0 3 or CA 6) precipitated in the slag skin were ph y s i c a l l y much thicker than they were in the slag bulk. In fact, the average content of aluminium in the slag skin i s much higher than in the nominal composition of the slag. It appears that the constant washing off of the secondary molten phase from the slag skin and the supply of fresh slag to the 69 region allows continuous p r e c i p i t a t i o n of the primary phase t i l l the passage between the c r y s t a l s becomes so constricted "that flow of l i q u i d is no longer possible. In t h i s scenario, the temperature of the slag skin-molten slag interface i s the liquidus temperature and the temperature of the c h i l l e d layer-"equilibrium" layer interface i s equal to the eutectic temperature (see f i g . 49). Using the values of thermal 1 8 0 0 J 1 7 0 0 1 6 0 0 1 5 0 0 J H 0 0 1 3 0 0 , 1 2 0 0 . • C o F , ' C o F 2 • A l 2 0 3 I C A g ) 2 0 3 0 4 0 5 0 m o s s • / . A l 2 0 3 Figure 49 - Phase diagram for CaF^-A^O^ conductivity given by Taylor and M i l l s et a l 6 1 a thickness of the "equilibrium" layer of 0.3 mm i s obtained, which f i t s in well with observations. 70 The high heat flux (nearly double that in the slag region) from the molten head of metal to the mould must s t i l l be explained. From the measurements of other workers, the molten steel i s generally at a temperature of 1600-1650°C at the ingot top surface periphery. For a 70:30 slag the liquidus i s about 1500°C. If conduction i s considered to be the major mode of heat transfer then the. difference in temperature at the slag skin hot face caused by the r i s e of molten metal i s not s u f f i c i e n t to explain the doubling of the heat flux. Radiation too i s not a sat i s f a c t o r y answer since the structure of the slag skin ensures that substantial scattering takes place, i . e . only an increase in e f f e c t i v e thermal conductivity may be expected. The only other reason could be a reduction of the a i r gap. On the face of i t , even th i s appears doubtful since there i s no reason why a couple of millimeters of molten metal head should do what 500 mm of slag head could not. However, i f i t i s supposed that the higher metal temperature w i l l melt the secondary phase and cause some dissolution of the primary c r y s t a l s , then the slag skin w i l l weaken s t r u c t u r a l l y without apparent change on subsequent r e - s o l i d i f i c a t i o n . This then, w i l l permit deformation of the slag skin and decrease the a i r gap. If the l i q u i d metal head i s excessively big, the skin may fracture and cause runout metal fi n s such as have been observed. 1'" 7 As the ingot s o l i d i f i e s and withdraws from the slag skin on contraction, the reduced pressure created draws out the secondary phase and th i s s o l i d i f i e s on the s o l i d ingot wall. 71 This explains why the t h i r d layer consisting of the secondary phase i s found in the slag skin. The spikes observed on t h i s surface mainly have e a s i l y oxidisable elements. It i s possible that corrosion products c o l l e c t near the metal meniscus since they are not readily soluble in the slag which may get trapped and stick to the slag skin thus creating the spikes. In the other slags examined, for example the 70:15:15 composition, there i s some melt back of the skin expected and th i s seems to be a reasonable explanation for the r i s e of the heat flux at the metal head l e v e l . Of course, some deformation of the slag skin i s also possible but t h i s does not change the situation q u a l i t a t i v e l y . This explanation of the v a r i a t i o n of the heat flux seems reasonable, but there i s s t i l l no way of estimating the heat flux. The thickness of the slag skin w i l l depend on the heat transfer c o e f f i c i e n t to the slag skin from the molten slag depths, t h i s in turn depending on the magneto-hydrodynamic flow f i e l d . The a i r gap w i l l depend on the thermo-physical properties of slag skins. Work in th i s area, i t seems, i s necessary. Some work of th i s nature has been reported by Medovar et a l 4 6 showing that 70:30 slags have a lower c o e f f i c i e n t of thermal expansion than the complex slags having MgO and S i 0 2 . A l 2 0 3 and CaO were found to increase the strength at the expense of reduced d u c t i l i t y while S i 0 2 raised the d u c t i l i t y of CaF 2-CaO-Al 20 3 system slags. 72 4.4 C r y s t a l l i s a t i o n In Slag Skins It i s useful, at t h i s juncture, to examine the c r y s t a l l i s a t i o n sequence according to the phase diagram. The diagram for the CaO-Al 20 3-CaF 2 system is given by M i l l s and Keene" 8 and i s shown in fig.50. The compositions of four common slags are marked out. The works of Z h m o i d i n " 9 ' 5 0 ' 5 1 ' 5 2 ' 5 3 and others show that there i s l i t t l e s o l i d s o l u b i l i t y and t h i s i s assumed in the following discussion. A 70:30 slag, in practice, contains up to 5% CaO due to hydrolysis of CaF 2 by moisture and evaporation of A1F 3. At high temperatures th i s slag i s in the 2-liquid region. The p l a i t points are unknown but from a knowledge of the slag skins and the work of Zhmoidin 5 0 they may be placed at points A and B as shown. The two l i q u i d s w i l l then have the composition C and D. If C has a higher liquidus the primary phase w i l l possibly be CA6 or A l 2 0 3 , both of which have been observed by M i t c h e l l 3 0 in the slag skin. This i s also supported by the work of Zhmoidin 5 0 who determined that CA6 would be the primary phase in t h i s system. A 70:15:15 slag l i e s at the edge of the 2 l i q u i d region though not in i t . The primary phase should be C 3A 3F with CaF 2 and CaO c r y s t a l l i s i n g l a t e r . However, small increases of CaO w i l l change the c r y s t a l l i s a t i o n sequence s i g n i f i c a n t l y . Therefore, i t i s l i k e l y that C^AyF w i l l come out instead of C 3A 3F. A 55:35:10 slag i s in the 2 l i q u i d region and following the arguments made e a r l i e r , the primary phase could be CA2 or C 3A 3F 73 mass V. CoO 70:30 70:15:15 55=35:10 33:33:33 F i g u r e 50 - P h a s e d i a g r a m f o r C a F 2 - A l ^ O g - C a O X JL m rr 74 depending on the t i e - l i n e s . Here, as for the 70:15:15 slag, small changes in the CaO content can change the p r e c i p i t a t i o n sequence. A 33:33:33 slag gave C^ i A 7 F with CaF 2 and CaO constituting the eutectic. M i t c h e l l 3 0 has found C 1 2A 7 and C 1 1 A 7 F in the slag skin of t h i s slag. Apparently c r y s t a l l i s a t i o n of the slag skin follows the phase diagram. Fractional c r y s t a l l i s a t i o n as spoken of by Kamensky et a l , 2 7 w i l l take place owing to the slag motion in the electroslag process as described above. This means that the chances of finding the secondary phases in the slag skin are low and one expects to find, b a s i c a l l y , only the primary phase in the bulk of the slag skin. There i s a great difference in the 70:30 slag owing to the c h a r a c t e r i s t i c s of c r y s t a l l i n e corundum (CA6 has a c l o s e l y similar structure to that of alpha A l 2 0 3 ) . This grows as a plate and so one finds CaF 2 in the skin in between the plates. As the plates thicken, one should find less of CaF 2 in the skin as may be seen in the slag skins of large ingots. 4.5 Slag Composition Effects On The Heat Flux It was suggested e a r l i e r that the metal temperature at the metal ingot top periphery i s between 1600-1650°C under normal melting circumstances. This temperature is above the melting point of the primary phases for a l l the slags mentioned in the preceding section with the exception of the 70:30 slag. 75 The implication i s that the slag skin w i l l be melted back in a l l cases except that of the 70:30 slag. The presence of l i g h t s t r i a t i o n s and the absence of "spikes" on the inside surface of a slag skin i s ind i c a t i v e of t h i s . The case for the 70:30 slag was dealt with in d e t a i l e a r l i e r . Continuing further, i f the slag skin i s thicker (as in large, ingot manufacture), a peak of heat flux may be observed at the molten metal head l e v e l in cases where the slag skin melts back and where i t does not melt, as for the 70:30 slag, the heat flux peak w i l l be d i f f u s e and may not even appear. This i s seen in the the results of Kondo et a l . 9 In s t a t i c mould conditions, the l a t t e r s i t u a t i o n favors a smoother surface. In the moving mould case, the system i s d i f f e r e n t . It has been observed that the skin of the 70:30 slag i s prone to fracture and thus defects, such as f i n s , are produced. Paton 5 8 indicates that the slag skin should be p l a s t i c to reduce the chances of fracture and suggests the addition of MgO and S i 0 2 to achieve t h i s . This " p l a s t i c i t y " i s probably due to the softening of glass present in the slag and the addition of S i 0 2 probably favors t h i s . It i s known that slags with low fl u o r i d e contents are operationally more e f f i c i e n t . 5 7 This i s due to their higher e l e c t r i c a l r e s i s t i v i t y and higher v i s c o s i t y which also means reduced slag flow v e l o c i t i e s . Kusamichi et a l , 8 report that there i s a discernable r a d i a l temperature gradient in slags of low fl u o r i d e content which i s absent in high fluoride slags. The net e f f e c t s are thicker slag skins and lower heat flux to 76 the mould for low fluoride slags. Further, from the discussion above, the slag skins of low f l u o r i d e (about 20-40%) slags have a higher melting point and so remelting of the slag skin should be l e s s . These skins should be stronger and so give better surfaces in moving mould operations. A 70:30 slag also w i l l give a thick skin but owing to the plate structure of corundum thi s skin fractures e a s i l y . The best skin, i t appears, must have a high melting point, must c r y s t a l l i s e with no d i r e c t i o n a l bias and should possess an appropriate degree of " p l a s t i c i t y " . 4.6 Factors Affecting Surface Quality From the discussion, i t appears that the mould cooling regime has l i t t l e to do with the surface q u a l i t y . One thing i s clear; the slow rate of r i s e of the ingot means that meniscus s o l i d i f i c a t i o n i s the only mechanism of formation of a bad rippled surface. Therefore, a good surface i s ensured by the presence of a l i q u i d metal head. Such a condition can be improved by reducing the resistance to heat flow in the l i q u i d metal and by increasing the resistance to heat flow to the mould. The former i s controlled by geometry and the power input to the system. The l a t t e r i s favoured by a thicker and/or stronger slag skin. A thicker skin can be generated by lower power input or by a thinner electrode, but t h i s i s counter productive in that the thermal resistance in the l i q u i d metal i s adversely 77 affected. The other method i s to select the correct slag. The use of pure compounds or eutectic compositions i s not recommended since in t h i s case remelting of the slag skin i s possible which causes a reduction of the thermal resistance. It i s necessary for a high melting phase to precipitate out and i t i s preferable that the proportion of the primary phase be as high as possible. Thus the tendency of the slag skin to dissolve and weaken i s combated. An increased slag depth at the same power input may also be b e n e f i c i a l but th i s i s a limited option in that a large increase of slag depth has a substantial penalty in power losses. 78 V. SUMMARY In conclusion, 1. Quasi steady-state three dimensional models may be used to calculate temperatures in mould sections to be used in electroslag casting. Temperature at the hot face may be changed by changing the spacing between channels but dr a s t i c changes can be made only by changing the mould material. For short production runs, temperature appears to be a v a l i d c r i t e r i o n , but for long runs design to minimise stresses may be necessary. 2. It appears that temperature d i s t r i b u t i o n s generated in t h i s way may be used to calculate thermal stresses, but v e r i f i c a t i o n i s necessary. 3. A proposal to explain the var i a t i o n of heat flux to the mould has been made. Work to determine i f the properties of slag skins can be used to estimate heat fluxes for arb i t r a r y shapes in a simple manner seems necessary. 4. The surface quality of a casting i s not influenced by mould design but by the slag composition, the power input and the geometrical relationship between the electrode, the slag depth and the casting shape. 79 REFERENCES 1. B.I.Medovar, V.L.Shevtsov, G.S.Marinsky, V.F.Demchenko, V.E.Machenko, Thermal Processes in Electroslag Flows , Naukova Dumka, Kiev, 1978. 2. A.Mitchell, R.M.Smailer, Int. Met. Rev., (5&6), 1979, p.231 . 3. W.Holzgruber, Stahl und Eisen, 88, (12), p.638 quoted in Ref. 4. 4. W.Duckworth, G.Hoyle, Electroslag Refining , Chapman and H a l l , London, 1969. 5. S.Joshi, Ph.D. thesis, University of B r i t i s h Columbia, 1 971 . 6. D.N.Pocklington, B.Patrick, Met. Trans., 10B, 1979, p.359-366. 7. Staff, Albany Metallurgy Research Centre, The Electr o s l a g Process , ed. R.H.Nafziger, U.S.Bureau of Mines, 1976, p.31-37. 8. T.Kusamichi, T . I s h i i , T.Onoye, K.Narita, Tetsu-to-Hagane, 12, 1980, p.1640-1649. 9. Y.Kondo, H.Kodama, E.Niyama, S.Hanada, H.Ishikawa, 101 ISIJ meeting, A p r i l 1981, lecture S256. 10. R.C.Sun, J.W.Pridgeon, Proc. Second Int. Symp. on ESR Technology, Pittsburgh, 1969. 11. L.F.Carvajal, G.E.Geiger, Met. Trans., 2, 1971, p.2087-2092. 12. J . F . E l l i o t t , M.Maulvault, Proc. Fourth Int. Symp. ESR Processes, Tokyo, 1973, p.69-80. 13. B.E.Paton, et a l . F i f t h Int. Symp. on ES and other Special Melting Technologies, Pittsburgh, 1974, p.323-344. 14. A.S.Ballantyne, Ph.D. thesis, University of B r i t i s h Columbia, 1978. 15. J.Kreyenberg, K.Schwerdtfeger, Arch. Eisen., 50, 1979, p. 1-6. 16. B.E.Paton, et a l , Special Electrometallurgy , Naukova Dumka, Kiev, 1972. 17. A.H.Dilawari, J.Szekely, Met. Trans., 11B, 1980, p.77-87. 80 18. M.Choudhary, J.Szekely, Met. Trans., 11B, 1980, p.439-453. 19. M.Choudhary, J.Szekely, Ironmaking and Steelmaking, 5, 1981, p.225-232. 20. R.S.Cremisio, E.D.Zak, Proc. Fourth Int. Symp. ESR Processes, Tokyo, 1973, p.137-148. 21. E.D.Zak, Proc. Vac. Met. Conf., Pittsburgh, 1977, p.315-331. 22. S.S.Kutaladze, V.M.Borishansky, A Concise Encyclopaedia of Heat Transfer , Pergamon Press, Oxford, 1966. 23. I.V.Samarasekera, J.K.Brimacombe, Int. Met. Rev., (6), 1978, p.286-300. 24. I.V.Samarasekera, J.K.Brimacombe, Sec. Process Tech. Conf., Chicago, 1981, p.2-21. 25. A.Mitchell, S.Bala, A.S.Ballantyne, G.Sidla, Metals Tech., Aug. 1981 , p306-312. 26. A.Mitchell, M.Etienne, Trans. AIME, 242, July 1968, p.1462-1463. 27. Y.M.Kamensky, Y.I.Afanasjev, B.N.Sakhotin, E.F.Vegoman, V.I.Yavoisky, L.A.Safronova, Reports of Int. Symp. on Special Electrometallurgy, Kiev, 1972. 28. B.Korousic, V.Osterc, Rad. Rundsch., 4, 1976, p.803-813. 29. M.Bell, A.Mitchell, JISI, 209, 1971, p.658-660. 30. A.Mitchell, unpublished. 31. K.Yamaguchi, M.Funazu, T.Ichihara, Proc. Fourth Int. Symp. ESR Processes, Tokyo, 1973, p.91-101. 32. P.O.Mellberg, Proc. Sixth Int. Vac. Met. Conf. on Special Melting, San Diego, 1979, p.535-542. 33. D.M.Longbottom, A.A.Greenfield, G.Hoyle, M.J.Rhydderch, Proc. Fourth Int. Symp. ESR Processes, Tokyo, 1973. 34. V.L.Shevtsov, B.I.Medovar, G.S.Marinski, V.P.Serdyukova, Henry Brutcher Trans. 9727, Refining Remelting , Naukova Dumka, Kiev, 1975, p.18-25. 35. A.Mitchell, G.Sidla, The Design, Construction and Operation of an ESC i n s t a l l a t i o n , Special Report to DREP/DSS, June 1980. 81 36. S.Patankar, Numerical Heat Transfer and F l u i d Flow , Hem i sphe r e(McGraw-Hill), 1980. 37. B.Carnahan, H.A.Luther, J.O.Wilkes, Applied Numerical Methods , J.Wiley and Sons, 1969. 38. F.Kreith, W.Z.Black, Basic Heat Transfer , Harper and Row, New York, 1980. 39. V.P.Simonov, et a l , BISITS trans. 12266, Izv. Vuz. Chern. Met., 1, 1974, p.175-178. 40. A.I.Chizhikov, V.L.Iokhimovich, G.P.Rachuk, L.I.Morozenskii, I.A.Mikhailova, S t a l ' in English, 11, 1968, p.921-924. 41. R.W.Young, V.Sunaryo, M.Cross, B r i t i s h Steel Corp. Report NO. Tech/722/1/75/A, July 1975. 42. R.Alberny, A.Leclerq, D.Amainy, M.Lahousse, Rev. de Met., 74, July-Aug 1976. 43. C.A.Sleicher, M.Tribus, Trans. ASME, 79, 1957, p.789-797. 44. E.M.Sparrow, R.SiegeL., Trans. ASME, 82, Aug 1960, p. 170. 45. V.J.Berry, App. S c i . Res., 4A, 1953, p.61-75. 46. B.I.Medovar, V.S.Zhakhovsky, A.G.Bogachenko, V.M.Matryn, Proc. Seventh ICVM, Tokyo, 1982, p.1244-1251. 47. G.Hoyle, Proc. Sixth Int. Vac. Conf. on Special Melting, San Diego, 1979, p.624-640. 48. K.C.Mills, B.J.Keene, Int. Met. Rev., (1), 1981, p.21-69. 49. G.S.Smirnov, A.K.Chatterjee, G.I.Zhmoidin, J. of Mat. S c i . , 8, 1973, p.1278-1282. 50. G.I.Zhmoidin, A.K.Chatterjee, Inorg. Mat., Proc. Akad. Nauk. USSR, 10, (10), 1974,p.1846-1851. 51. A.K.Chatterjee, G.I.Zhmoidin, INorg. Mat., Proc. Akad. Nauk. USSR, 8, (5), 1972, p.886-892. 52. G.I.Zhmoidin, G.S.Smirnov, V.N.Glotky, Inorg. Mat., Proc. Akad. Nauk. USSR, 13, (3), 1977, p.552-553. 53. G.I.Zhmoidin, Metals, Proc. Akad. Nauk. USSR, 6, 1969, p.9-15. 54. V.M.Borishansky, I.I.Paleev, A.A.Andreevsky, B.S.Fokin, 82 M.E.Lavzentiev, K.P.Malyus-Malitsky, V.N.Fromzel, G.P.Danilova, Int. J . Heat Mass Trans., (5), 1973, p. 1073-1085. 55. E.N.Sieder, G.N.Tate, Ind. Eng. Chem., 28, (12), 1936, p.1429. 56. B.S.Petukhov, Adv. in Heat Transfer, 6, 1970, p.504-564. 57. J.D.W.Rawson, G.Jeszensky, A.W.Bryant, Proc. Sixth Int. Vac. Met. Conf. on Special Melting, San Diego, 1979, p.848-863. 58. B.E.Paton, B.I.Medovar, V.L.Artamonov, A.G.Bogachenko, Y.P.Shtanko, Henry Brutcher Trans. 9726, Refining Remelting , Naukova Dumka, Kiev, 1972, p.49-53. 59. K.C.Mills, J.S.Powell, J.W.Bryant, B.J.Keene, Can. Met. Quart., 20, (1), 1981, p.93-99. 60. A.A.Sharapov, C.E.Volkov, C.V.Lebedev, Teoriya Metal., 3, 1957, p.131-136. 61. R.Taylor, K.C.Mills, Arch. Eisen., 53, (2), 1982, p.55-63. 62. R.N.Noyes, Trans. ASME, 83, 1961, p.96-97. 83 APPENDIX A - SLAG NOTATION i . 70:30 = 70% CaF 2 , 30% A l 2 0 3 i i . 70:15:15 = 70% CaF 2, 15% A l 2 0 3 , 15% CaO. i i i . 50:30:20 = 50% CaF 2, 30% A l 2 0 3 , 20% CaO. i v . 33:33:33 = 1/3 CaF 2, 1/3 A l 2 0 3 , 1/3 CaO. 84 APPENDIX B ~ MISC. DATA AND CORRELATIONS i . Water flow measurements through Al section 30 l i t r e s / m i n . through Cu section 30 l i t r e s / m i n . through heat flux sensor 14 l i t r e s / m i n . These measurements were made by noting the time needed to f i l l a 100 l i t r e drum. i i . Some physical properties Material Thermal Conductivity (W/m.K) at 273K at 400K at 600K Aluminium 236 240 232 Copper 401 392 383 Water at sat. Density Absolute Prandtl Pr. and at ViscosityNo. kg/m3 X10 6 m2/s 273 K 999.3 1794 13.7 313 K 999.2 658 4.3 353 K 971.8 352 2.25 373 K 958.4 278 1.75 i i i . Sieder-Tate c o r r e l a t i o n Nu = .027 Re 8 Pr 3 3 (nb /ns ) •1 • where Nu = Nusselt number Re = Reynolds number Pr = Prandtl number Mi = v i s c o s i t y of f l u i d at bulk temperature Mj = v i s c o s i t y of f l u i d at surface of tube 85 A l l properties except v i s c o s i t y taken at average bulk temperature. i v . Petukhov co r r e l a t i o n Nu = U /8) Re Pr K , U ) + K 2(Pr) U/8) -5 (Pr ' 6 7 - 1) where £ = (1.82 Log Re - 1.64)" 2 K, = 1 + 3.4 £ K 2 = 11.7 + 1.8 P r " • 3 3 A l l properties taken at the mean f i l m temperature 86 APPENDIX C ~ FORTRAN PROGRAM TO CALCULATE TEMPERATURES IN CHANNEL COOLED MOULDS C c C THIS IS A PROGRAMME TO CALCULATE TEMPERATURES IN C CHANNEL COOLED MOULDS USED IN ELECTROSLAG CASTING C THE METHOD IS TO USE MATRIX METHODS TO CALCULATE C THE TEMPERATURES FOR EACH SLICE AND ITERATE FROM C SLICE TO SLICE TILL CONVERGENCE IS REACHED. C THE REQUIRED INPUT ON UNIT 5 IS AS GIVEN BELOW. C THE INITIAL TEMP. ESTIMATE MUST BE GIVEN ON UNIT 4 C AND THE OUTPUT IS GIVEN ON UNIT 7. Q ******************************* C M = NUMBER OF 'X' DIVISIONS C N = NUMBER OF Y NODES C 0 = NUMBER OF Z NODES C L = MAX. NO. OF ITERATIONS ALLOWED C TOLER = TOLERANCE C TW = INLET TEMPERATURE OF COOLANT C H = HEAT TRANSFER COEFFICIENT TO THE COOLANT C DELX = DISTANCE BETWEEN 2 X NODES C DELZ = DISTANCE BETWEEN 2 Z NODES C RH = RADIUS OF ROUND CHANNEL IN TERMS OF DELX UNITS C CH = POSITION OF CHANNEL IN TERMS OF DELX UNITS C COTC = COEFFICIENT OF THERMAL CONDUCTIVITY AT 0 C C VELA = MASS FLOW VELOCITY OF COOLANT C CEE = GRADIENT OF THERMAL CONDUCTIVITY WITH TEMP. C QA(I) = HEAT FLUX INTO THE HOT FACE C I N I T = 0 I F THERE IS NO GUESS SOLUTION, OTHERWISE ANY NO. C IT IS SUGGESTED THAT ALL UNITS BE FED IN IN CGS UNITS C TO AVERT FORMAT PROBLEMS. C ********************************************************** c C [TE] {S} = {B} IS THE MATRIX EQUATION TO BE SOLVED C AT EACH SLICE. X,Y,Z ARE SPECIAL VARIABLES. C T — TEMPERATURE C TT -- REFERENCE TEMPERATURE FIELD C I PERM, TS — MEMORY REQUIRED BY SOLUTION ROUTINE Q ****************************************************** C THE ROUTINE FSLE IS A ROUTINE IN THE UBC LIBRARY AND C USES GAUSSIAN ELIMINATION TO SOLVE SIMULTANEOUS C EQUATIONS. C ************************************************************** C THE ROUTINES MATDAT, BOTDAT AND TOPDAT SERVE TO CALCULATE C THE MATRIX TERMS FOR THE MIDDLE SLICES, THE BOTTOM SLICE C AND THE TOP SLICE RESPECTIVELY C c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * REAL*4 TT,B,T,DET DIMENSION TT(20,20,30),TS(400,400),IPERM(800),S(400),X(20) 87 1,Y(20),HT(30) COMMON T(20, 20,30),TE(400,400),B(400),QA(30),Z(20 , 20) , 1TH(30),DELT C C READ IN MATRIX DATA C INTEGER 0 READ(5,1) M,N,0,L,TOLER 1 FORMAT(2X,3I3,I4,F3.1) READ(5,2) TW,H 2 FORMAT(2X,F4.1,F7.3) READ(5,3) DELX,RH,CH,DEL Z 3 FORMAT(2X,F6.4,2F4.1,F6.4) READ(5,4) COTC,VELA,CEE 4 FORMAT(2X,F5.4,F7.3,F8.6) READ(5,6) (QA(I),1=1,0) 6 FORMAT(1X,10F6.2) READ(5,7) INIT 7 FORMAT(1X,I 2) IF (INIT .EQ. 0) GO TO 9 CONTINUE READ(4,8) (((T(I,J,K),J=1,N),I=1,M),K=1,0) 8 FORMAT(11F6.2) C ********************************************* C THE LINES BELOW TO BE USED TO VARY H ALONG THE HEIGHT OF C THE MOULD. Q *********************************************************** C H T ( 1 ) = H - 0 . 1 8 C H T ( 2 ) = H - 0 . 1 6 C HT(3) = H- 0.16 C HT(4) = H-0.14 C HT(5) = HT(4) C H T ( 6 ) = H - 0 . 1 0 C HT(7) = H - 0.08 C HT(8) = H C HT(9) = H+.02 C HT(10) = H+.04 C HT(11) = H +.02 C HT(12) = H- .06 C HT(13) = H - 0.08 C HT(14) = H-0.10 C HT(15) = H-0.12 C HT(16) = H - 0.14 C HT(17) = H- 0.16 C HT(18) = HT(17) Q *********************************************************** MN=M*N C C GIVE X, Y COORDINATES TO THE NODES C DO 10 1=1,M DO 10 J=1,N DO 10 K=1,0 10 T(I,J,K) = 0.0 88 SUM =-1.0 DO 11 1 = 1 ,M SUM=SUM+1. 11 X(I) = SUM SCUM =-1.0 DO 12 J=1 ,N SCUM=SCUM+1. 12 Y(J) = SCUM C C CODE DIFFERENT TYPES OF NODES C DO 100 1 = 1 ,M DO 100 J=1 ,N C C THE FOLLOWING IS BASED ON COORDINATE GEOMETRY C AB =ABS(Y(J)-CH) A = (X(I)**2 + (Y(J) - CH)**2) IF ((RH**2) .LE. A) GO TO 20 C C THE NODE IS IN THE COOLING CHANNEL C Z(I,J) = 1 GO TO 100 20 IF (X(I) .GE. (RH+1.)) GO TO 50 IF (ABS(Y(J)-CH) .GE. (RH+1.)) GO TO 50 D = ABS(SQRT(ABS(RH**2-X(I)**2))-AB) C = ABS(X(I) - SQRT(ABS(RH**2 - (AB)**2))) IF (D .GE. 1.) GO TO 28 IF (X(I).GT. 0.) GO T022 IF (Y(J) .GT.CH) GO TO 21 C C CODES 2 TO 9 REFER TO NODES BOUNDING THE COOLING CHANNEL C Z(I,J) = 2 GO TO 100 21 Z(I,J) =3 GO TO 100 22 IF (C.GE.1.) GO TO 26 IF (Y(J)-CH) 23,24,25 23 Z(I,J) = 4 GO TO 100 24 Z(I,J) = 5 GO TO 100 25 Z(I,J) = 6 GO TO 100 26 IF (Y(J).GT.CH) GO TO 27 Z(I,J) =7 GO TO 100 27 Z(I,J) = 8 GO TO 100 28 IF (C.GE.1.) GO TO 50 Z(I,J) = 9 GO TO 100 89 50 IF (X(I).LT.(M-1)) GO TO 65 IF (Y(J).GT.O.) GO TO 55 C " C CODES 10 TO 12 REFER TO THE LEFT SIDE OF THE SECTION C INCLUDING THE CORNER NODES BUT EXCLUDING THE NODES C INSIDE THE COOLING CHANNEL AND THOSE BOUNDING IT. C Z(I,J) = 10 GO TO 100 55 IF (Y(J).LT.(N-1)) GO TO 60 Z(I,J) = 11 GO TO 100 60 Z(I,J) = 12 GO TO 100 65 IF (X(I).GT.O.) GO TO 80 IF (Y(J).GT.O.) GO TO 70 C C CODES 13 TO 15 REFER TO THE NODES ON THE RIGHT SIDE OF C THE SECTION INCLUDING THE CORNER NODES C Z(I,J) = 13 GO TO 100 70 IF (Y(J).LT.(N-1)) GO TO 75 Z(I,J) =14 GO TO 100 75 Z(I,J) =15 GO TO 100 80 IF (Y(J).GT.O.) GO TO 85 C C CODE 16 REFERS TO THE BOUNDARY OPPOSITE THE HOT FACE C Z(I,J) =16 GO TO 100 85 IF (Y(J).LT.(N-1)) GO TO 90 C C CODE 17 REFERS TO THE HOT FACE BOUNDARY C Z(I,J) =17 GO TO 100 C C CODE 18 REFERS TO INTERIOR NODES C 90 Z(I,J) =18 100 CONTINUE C C SET THE TEMP. OF THE CHANNEL TO THAT OF THE COOLANT C AT THE INLET C DO 101 IC=1,0 101 TH(IC) = TW C C SET ASSUMED TEMP. FIELD = REFERENCE C DO 105 K=1,0 90 DO 105 J=1,N DO 105 1 = 1 ,M 105 TT(l,J,K) = T(I,J,K) C C ITERATE L TIMES C DO 310 LL=1,L C C RESET INCREASE OF TEMP. OF THE COOLANT TO ZERO C FOR THE NEXT SLICE C DELT =0.0 Q *************************************** C TO BE USED ONLY IF H IS TO BE VARIED C C DO 130 K=1,0 C H = HT(K) c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C BOTTOM SLICE C IF (K.EQ.1) GO TO 110 C C TOP SLICE C IF (K.EQ.O) GO TO 120 CALL MIDDAT(M,N,O,K,COTC,DELX,DELZ,VELA,H,CEE) CALL FSLE(MN,400,TE., 1 , 400 , B, S , I PERM, 400 , TS , DET, JEXP) C C RE-WRITE SOLUTION IN TERMS OF X,Y,Z COORDINATES C DO 106 11=1,MN I = (H-O/N+1 J= II-N*(I-1) 106 T(I,J,K) = S(II) GO TO 130 110 CALL BOTDAT(M,N,0,K,COTC,DELX,DELZ,VELA,H,CEE) CALL FSLE(MN,400,TE, 1 ,400,B,S,I PERM,400,TS,DET,JEXP) DO 115 11=1,MN I=(II-1)/N+1 J= II-N*(I-1) 115 T(I,J,K) = S(II) GO TO 130 120 CALL TOPDAT(M,N,0,K,COTC,DELX,DELZ,VELA,H,CEE) CALL FSLE(MN,400,TE,1,400,B,S,I PERM,400,TS,DET,JEXP) DO 125 11=1,MN I = (II-1)/N+1 J = II-N*(I-1) 125 T(I,J,K) = S(II) 130 CONTINUE C C COMPARE REFERENCE AND CALCULATED FIELDS C 91 DO 200 1=1,M DO 200 J=1,N DO 200 K=1,0 IF (ABS(TT(I,J fK)-T(l,J,K)) .GT. TOLER) GO TO 300 200 CONTINUE C C WRITE OUT RESULTS. NOTE IF N > 20 THEN PRINTOUT WILL C LOOK FUNNY C DO 218 K=1,O DO 216 1 = 1 ,M WRITE(7,210) (T(I,J,K),J=1,N) 210 FORMAT(20F7.2) 216 CONTINUE WRITE(7,217) 217 FORMAT(/,/) 218 CONTINUE STOP C C RESET REFERENCE FIELDS C 300 DO 310 K=1,0 DO 310 J=1,N DO 310 1=1,M 310 TT(I,J,K) = T(I,J,K) C C NO. OF ITERARIONS WERE INSUFFICIENT. WRITE OUT FAILURE C AND TEMPERATURES. C WRITE(7,320) TOLER,L 320 FORMAT(1H ,'THE TEMPERATURES DO NOT CONVERGE TO WITHIN' 1' DEGREES IN',14,' ITERATION STEPS.') DO 328 K=1,0 DO 326 1=1,M WRITE(7,321) (T(I,J,K),J=1,N) 321 FORMAT(20F7.2) 326 CONTINUE WRITE(7,327) 327 FORMAT(/,/) 328 CONTINUE STOP END SUBROUTINE MIDDAT(M,N,O,K,COTC,DELX,DELZ,VELA,H,CEE) C C THIS ROUTINE CALCULATES MATRIX TERMS FOR THE MIDDLE SLICES C COMMON T(20,20,30),TT(400,400),B(400),QA(30),Z(20,20), 1TH(30),DELT INTEGER CODE,0 MN=M*N C C THE FOLLOWING CALCULATIONS ARE EQUIVALENT TO THE VARIOUS C COEFFICIENTS. THEY ARE MADE HERE TO AVOID REATING AGAIN C 92 H12 = H*DELX*DELZ*.5 H1 = H12*2. H2 = H12*4. H12TH = H12*TH(K) H1TH= H12TH*2. H2TH = H12TH*4. QA12 = QA(K)*DELX*DELZ*.5 QA1 = QA12*2. C C INITIALISE MATRICES TO ZERO C CALL GSET(TT,400,400,400,0.) CALL GSET(B,400,1 ,400,0. ) THEAT =0.0 DDH = DELZ * DELX*H C C DO FOR EACH NODE DEPENDING ON CODE C DO 100 J=1,N DO 100 1=1,M C C SET THERMAL CONDUCTIVITY ACCORDING TO TEMP. OF NODE C AS IN THE PREVIOUS ITERATION C CT = CEE*T(I,J,K) + COTC COT12 = CT*.5*DELZ COT1 = COT12*2. COT2 = COT1*2. COT3 = COT1*3. COT4 =COT1*4. C14 = CT*(DELX**2)*.25/DELZ C12 = C14*2. C = C14*4. CODE = Z(I,J) KK = N*(I-1)+J GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),CODE 1 TT(KK,KK) = 1. B(KK) = TH(K) GO TO 100 2 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT1+ H12 + C12) TT(KK,KK+N) = COT12 B(KK) = -(H12TH+ C14*(T(I,J,K+1) + T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT*.5 GO TO 100 3 TT(KK,KK) = -(COT1+ H12+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(H12TH+ C14*(T(I,J,K+1) + T(I,J,K-1))) HEAT = DDH*( T( I , J , K ) - T H ( K ) ) THEAT = THEAT + HEAT * .5 GO TO 100 4 TT(KK,KK-1) = COT1 93 TT(KK,KK) = -(C0T2+ H2+ C*2.) TT(KK,KK+N) = C0T1 B(KK) = -(H2TH+ C*(T(I,J,K+1) + T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT+ HEAT *2. GO TO 100 5 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ H1+ C) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(H1TH + C12*(T(I,J,K+1)+T(l,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 1. GO TO 100 6 TT(KK,KK) = -(COT2+ H2+ C*2.) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT1 B(KK) = -(H2TH+ C*(T(I,J,K+1)+T(I,J,K~1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT*2. GO TO 100 7 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ H1+ C) TT(KK,KK-N) = COT12 TT(KK,KK+N) = COT12 B(KK) = ~(H1TH+C12*(T(I,J,K+1)+T(l,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT*1. GO TO 100 8 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT2+ H1 + C) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT12 B(KK) = -(H1TH + C12*(T(I,J,K+1) + T( I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT+ HEAT*1. GO TO 100 9 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ H1+ C) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(H1TH+ C12*(T(I,J,K+1) + T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 1. GO TO 100 10 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT1+ C12) TT(KK,KK+1) = COT12 B(KK) = -(C14*(T(I,J,K+1) + T(I,J,K-1))) GO TO 100 1 1 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT1+ C12) TT(KK,KK-1) = COT1 2 B(KK) = -(QA12+ C14*(T(I,J,K+1) + T(I,J,K-1))) 94 GO TO 100 12 TT(KK,KK-N) = COT1 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ C) TT(KK,KK+1) = COT1 2 B(KK) = -(C12*(T(I,J,K+1) + T(I,J,K-1))) GO TO 100 13 TT(KK,KK) = -(COT1+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(C14*(T(I,J,K+1) + T(I,J,K-1))) GO TO 100 14 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT1+ C12) TT(KK,KK+N) = COT12 B(KK) = -(QA12+ C14*(T(l,J,K+1) + T(I,J,K-1))) GO TO 100 15 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ C) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(C12*(T(I fJ,K+1) + T ( I,J,K-1))) GO TO 100 16 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT2+ C) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT12 B(KK) = -(C12*(T(I,J,K+1) + T ( I,J,K-1))) GO TO 100 17 TT(KK,KK-N) = COT12 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ C) TT(KK,KK+N) = COT12 B(KK) = -(QA1+ C12*(T(I,J,K+1) + T ( I,J,K-1))) GO TO 100 18 TT(KK,KK-N) = COT1 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT4+ C*2.) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT1 B(KK) = -(C*(T(I,J,K+1) + T(I,J,K-1))) 100 CONTINUE C C CALCULATE INCREASE OF TEMPERATURE AND ADD TO TEMP. C OF COOLANT C DELT = DELT + THEAT/VELA IF (DELT .LT. 0.001) GO TO 120 K01=K+1 DO 110 KO=K01,0 110 TH(KO) = TH(1) + DELT 120 CONTINUE RETURN END 95 SUBROUTINE BOTDAT(M,N,0,K,COTC,DELX,DEL Z,VELA,H,CEE) C C THIS ROUTINE SETS MATRIX TERMS FOR THE BOTTOM SLICE C ALL COMMENTS GIVEN IN THE ROUTINE MIDDAT APPLY HERE TOO C COMMON T(20,20,30),TT(400,400),B(400),QA(30),Z(20,20), 1TH(30),DELT INTEGER CODE,0 MN=M*N H12 = H*DELX*DELZ*.25 H1 = H12*2. H2 = H12*4. H12TH = H12*TH(K) H1TH= H12TH*2. H2TH = H12TH*4. QA12 = QA(K)*DELX*DELZ*.25 QA1 = QA12*2. CALL GSET(TT,400,400,400,0.) CALL GSET(B,400,1,400,0.) THEAT = 0.0 DDH = DELZ*DELX*H DO 100 J=1,N DO 100 1 = 1 ,M CT= CEE * T(I,J,K) +COTC COT12 = CT*.25*DELZ COT1 = COT12*2. COT2 = COT1*2. COT3 = COT1*3. COT4 =COT1*4. C14 = CT*(DELX**2)*.25/DELZ C12 = C14*2. C = C14*4. CODE = Z(I,J) KK = N*(I-1)+J GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),CODE 1 TT(KK,KK) = 1. B(KK) = TH(K) GO TO 100 2 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT1+ H12 + C14) TT(KK,KK+N) = COT12 B(KK) = -(H12TH+ C14*(T(I,J,K+1) + 0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT *.25 GO TO 100 3 TT(KK,KK) = -(COT1+ H12+ C14) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(H12TH+ C14*(T(I,J,K+1) + 0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.25 GO TO 100 4 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ H2+ C) 96 TT(KK,KK+N) = C0T1 B(KK) = -(H2TH+ C*(T(I,J,K+1) + 0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 1. GO TO 100 5 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ H1+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(H1TH + C12*(T(I,J,K+1)+0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.5 GO TO 100 6 TT(KK,KK) = -(COT2+ H2+ C) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT1 B(KK) = -(H2TH+ C*(T(I,J,K+1)+0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 1. GO TO 100 7 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ H1+ C12) TT(KK,KK-N) = COT12 TT(KK,KK+N) = COT 12 B(KK) = -(H1TH+C12*(T(I,J,K+1)+0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT *.5 GO TO 100 8 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT2+ H1 + C12) TT(KK,KK+1) = COT1 TT(KK,KK+N) - COT12 B(KK) = -(H1TH + C12*(T(I,J,K+1) + 0.)) HEAT = DDH*(T(I,J,K)-TH(K)) THEAT = THEAT + HEAT * 0.5 GO TO 100 9 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ H1+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(H1TH+ C12*(T(I,J,K+1) + 0.)) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT *0.5 GO TO 100 10 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT1+ C14) TT(KK,KK+1) = COT12 B(KK) = -(C14*(T(I,J,K+1) + 0.)) GO TO 100 1 1 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT1+ C14) TT(KK,KK-1) = COT12 B(KK) = -(QA12+ C14*(T(I,J,K+1) + 0.)) GO TO 100 97 12 TT(KK,KK-N) = C0T1 TT(KK,KK-1) = C0T12 TT(KK,KK) = -(C0T2+ C12) TT(KK,KK+1) = C0T12 B(KK) = -(C12*(T(I,J,K+1) + 0.)) GO TO 100 13 TT(KK,KK) = -(COT1+ C14) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(C14*(T(I,J,K+1) + 0.)) GO TO 100 14 TT(KK,KK-1) = COT1 2 TT(KK,KK) = -(COT1+ C14) TT(KK,KK+N) = COT12 B(KK) = -(QA12+ C14*(T(I,J,K+1) + 0.)) GO TO 100 15 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(C12*(T(I,J,K+1) + 0.)) GO TO 100 16 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT12 B(KK) = -(C12*(T(I,J,K+1) + 0.)) GO TO 100 17 TT(KK,KK-N) = COT12 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+N) = COT12 B(KK) = -(QA1+ C12*(T(I,J,K+1) + 0.)) GO TO 100 18 TT(KK,KK-N) = COT1 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT4+ C) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT1 B(KK) = -(C*(T(I,J,K+1) + 0.)) 100 CONTINUE DELT = DELT + THEAT/VELA IF (DELT .LT. 0.001) GO TO 120 K01=K+1 DO 110 KO=K01,0 110 TH(KO) = TH(1) + DELT 120 CONTINUE RETURN END SUBROUTINE TOPDAT(M,N,O,K,COTC,DELX,DEL Z,VELA,H,CEE) C C THIS ROUTINE CALCULATES MATRIX TERMS FOR THE TOP SLICE C ALL OTHER COMMENTS GIVEN IN THE ROUTINE MIDDAT APPLY C 98 COMMON T(20,20,30),TT(400,400),B(400),QA(30),Z(20,20), 1TH(30),DELT INTEGER CODE,0 MN=M*N H12 = H*DELX*DELZ*.25 H1 = H12*2. H2 = H12*4. H12TH = H12*TH(K) H1TH= H12TH*2. H2TH = H12TH*4. QA12 = QA(K)*DELX*DELZ*.25 QA1 = QA12*2. CALL GSET(TT,400,400,400,0.) CALL GSET(B,400,1,400,0.) THEAT = 0.0 DDH = 0.0 DO 100 J=1 ,N DO 100 1 = 1 ,M CT = CEE * T(I,J,K) +COTC COT12 = CT*.25*DELZ COT1 = COT12*2. COT2 = COT1*2. COT3 = COT1*3. COT4 =COT1*4. C14 = CT*(DELX**2)*.25/DELZ C12 = C14*2. C = C14*4. CODE = Z(I,J) KK = N*(I-1)+J GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),CODE 1 TT(KK,KK) = 1. B(KK) = TH(K) GO TO 100 2 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT1+ H12 + C14) TT(KK,KK+N) =.COT12 B(KK) = -(H12TH+ C14*(0. + T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.25 GO TO 100 3 TT(KK,KK) = -(COT1+ H12+ C14) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(H12TH+ C14*(0. + T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT *0.25 GO TO 100 4 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ H2+ C) TT(KK,KK+N) = COT1 B(KK) = -(H2TH+ C*(0. + T(I,J,K~1))) HEAT = DDH*(T(I,J,K) -TH(K)) THEAT =THEAT + HEAT * 1 . GO TO 100 99 5 TT(KK,KK-1) = C0T12 TT(KK,KK) = -(C0T2+ H1+ C12) TT(KK,KK+1) = C0T12 TT(KK,KK+N) = C0T1 B(KK) = -(H1TH + C12*(0.+T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.5 GO TO 100 6 TT(KK,KK) = -(COT2+ H2+ C) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT1 B(KK) = -(H2TH+ C*(0.+T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 1 . GO TO 100 7 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ H1+ C12) TT(KK,KK-N) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(H1TH+C12*(0.+T(I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.5 GO TO 100 8 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT2+ H1 + C12) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT12 B(KK) = -(H1TH + C12*(0. + T(I,J,K-1))) HEAT = DDH * (T(l,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.5 GO TO 100 9 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ H1+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(H1TH+ C12*(0. + T ( I,J,K-1))) HEAT = DDH*(T(I,J,K) - TH(K)) THEAT = THEAT + HEAT * 0.5 GO TO 100 10 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT1+ C14) TT(KK,KK+1) = COT12 B(KK) = -(C14*(0. + T ( I,J,K-1))) GO TO 100 11 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT1+ C14) TT(KK,KK-1) = COT12 B(KK) = -(QA12+ C14*(0. + T(I,J,K-1))) GO TO 100 12 TT(KK,KK-N) = COT1 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+1) = COT12 B(KK) = -(C12*(0. + T(I,J,K-1))) 100 GO TO 100 13 TT(KK,KK) = -(COT1+ C14) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT12 B(KK) = -(C14*(0. + T(I,J,K-1))) GO TO 100 14 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT1+ C14) TT(KK,KK+N) = COT12 B(KK) = -(QA12+ C14*(0. + T(I,J,K-1))) GO TO 100 15 TT(KK,KK-1) = COT12 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+1) = COT12 TT(KK,KK+N) = COT1 B(KK) = -(C12*(0. + T(I,J,K-1))) GO TO 100 16 TT(KK,KK-N) = COT12 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT12 B(KK) = -(C12*(0. + T(I,J,K-1))) GO TO 100 17 TT(KK,KK-N) = COT12 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT2+ C12) TT(KK,KK+N) = COT12 B(KK) = -(QA1+ C12*(0. + T(I,J fK-1))) GO TO 100 18 TT(KK,KK-N) = COT1 TT(KK,KK-1) = COT1 TT(KK,KK) = -(COT4+ C) TT(KK,KK+1) = COT1 TT(KK,KK+N) = COT1 B(KK) = -(C*(0. + T(I,J,K-1))) 100 CONTINUE DELT = THEAT/VELA TH(K) = TH(K) RETURN END 101 APPENDIX D - FORTRAN PROGRAM REFERRED TO IN SECTION 3.3 Q ********************************************** IMPLICIT REAL*8(A-H,0-Z) INTEGER BCC C THIS IS A PROGRAM USING TRIANGULAR OR RECTANGULAR FINITE C ELEMENTS TO SOLVE STEADY STATE HEAT TRANSFER PROBLEMS. C SPECIFICALLY IT IS BEING USED TO MODEL SECTIONS IN AN C ELECTROSLAG CASTING MOULD. ALL VARIABLES BEGINNING WITH C THE LETTERS A~H AND 0~Z ARE DECLARED DOUBLE PRECISION. C C DEFINE VARIABLES C C = CONDUCTANCE MATRIX C T = TEMPERATURE MATRIX C B = [C]*[T] MATRIX C ICO = ELEMENT CORNER"NODE NUMBERS C X,Y = GLOBAL CO-ORDINATES C D = CONDUCTIVITY MATRIX C H = HEAT TRANSFER CO-EFFICIENT( C ITYPE = TYPE OF ELEMENT C IMAT = MATERIAL CODE C BCC = BOUNDARY CODE(CONVECTION) C BCT = BOUNDARY CODE (TEMPERATURE) C IHF = CODE TO CALCULATE HEAT FLUX C XX,YY,A,S ELEMENT CO-ORDINATES AND MATRICES C C C DIMENSION C(200,200),T(200),B(200),X(200),Y(200),ITYPE(200) 1BCT(200),IHF(200),IMAT(200),IPERM(400),TT(200,200),DA(2,2), 2XX(8),BCC(200),YY(8),NO(8),ICO(200,8),S4(4,4),S3(3,3),A4(4) 3S8(8,8),A8(8),A3(3),DC(2,2),RZ(200) COMMON /ZD/ DET,JEXP COMMON /DSLMP$/ NITER C C CALL SUBROUTINE TO READ IN NECESSARY DATA FOR PROBLEM C CALL LAYOUT(X,Y,ICO,2 0 0,1TYPE,NE,NNODES,BCT,BCC,H,TW,COTC, 1IHF,COTA,IMAT) C C INITIALISE ARRAYS C CALL DGSET(DC,2,2,2,0.D0) DC(1,1) = COTC DC(2,2) = DC(1,1) CALL DGSET(DA,2,2,2,0.D0) DA(1,1) = COTA DA(2,2) = DA(1,1) DO 3 1=1,200 B(I) = 0.D0 DO 3 J=1,200 102 3 C(I,J) = 0.D0 C C BEGIN DO-LOOP TO GO THROUGH ALL ELEMENTS C DO 8 IEL=1,NE NN = ITYPE(IEL) BC = BCC(IEL) DO 4 I=1,NN NO(I) = IC0(IEL,I) XX(I) = X(NO(D) 4 YY(I) = Y(NO(l)) C C CALL SUBROUTINE TO CALCULATE I-TH CONDUCTANCE MATRIX C IM = IMAT(IEL) IF (IM .EQ. 0) GO TO 20 CALL COND(IEL,NN,XX,YY,S 3,S 4,S 8,DC,A3,A4,A8,H,TW,BC) GO TO 21 20 CALL C0ND(IEL,NN,XX,YY,S3,S4,S8,DA,A3,A4,A8,H,TW,BC) 21 CONTINUE C C CALL SUBROUTINE TO INSERT ELEMENT MATRIX INTO STRUCTURE C MATRIX C IF (NN .EQ. 3) GO TO 5 IF (NN .EQ. 8) GO TO 6 CALL SETUP(NO,NN,C,200,S4,NN,B,A4) GO TO 7 5 CALL SETUP(NO,NN,C,200,S3,NN,B,A3) GO TO 7 6 CALL SETUP(NO,NN,C,200,S8,NN,B,A8) 7 CONTINUE 8 CONTINUE C C CALL SUBROUTINE TO ADJUST STRUCTURE MATRIX FOR BOUNDARY C CONDITIONS C DO 50 1=1,NNODES IF (BCT(I)) 50,40,40 40 CALL BOUND(l,BCT,C,B,200,NNODES) 50 CONTINUE WRITE(7,55) ((C(I,J),1=1,NNODES),J=1,NNODES) 55 FORMAT(1X,2D30.22) WRITE(7,56) (B(I),1=1,NNODES) 56 FORMAT(/,/,/,IX,(2D30.12)) DO 57 1=1,NNODES B(I) = B(I)*1.D45 DO 57 J=1,NNODES 57 C(I,J) = C(I,J) * 1.D45 C C CALL MATRIX SOLVER FROM UBC LIBRARY C DEPS = 1.D-16 CALL DSLIMP(C,TT,B,T,RZ,I PERM,NNODES,200,DEPS,1,14) 103 WRITE(6,25) DET,JEXP 25 FORMAT(1X,' DETERM. = ',G20.12,'* 10**' ,I 5,/) ' WRITE(6,26) ((IJK,T(IJK)),IJK=1,NNODES) 26 F0RMAT(1X,4('T(',13,') = ',F6.2,3X)) 27 CONTINUE C C CALCULATE HEAT FLUX AT MOLD INSIDE SURFACE AND PRINT C DO 100 IEL=1,NE IK = IHF(IEL) IF (IK .EQ. 0) GO TO 100 NN = ITYPE(IEL) DO 90 1=1,NN NO(I) = IC0(IEL,I) XX(I) = X(NO(l)) YY(I) = Y(NO(l)) T(I) = T(NO(l)) 90 CONTINUE IM = IMAT(IEL) IF (IM .EQ. 0) GO TO 95 CALL HFLUX(IEL,NN,XX,YY,T,DC,NO) GO TO 96 95 CALL HFLUX(lEL,NN,XX,YY,T,DA,NO) 96 CONTINUE 100 CONTINUE STOP END Q ************************************** SUBROUTINE COND(l,NN,X,Y,S3,S4,S8,D,A3,A4,A8,H,TW,BC) IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(1),Y(1),S3(3,3),S4(4,4),D(2,2) 1,S8(8,8),A8(8),A3(3),A4(4) IF (NN.EQ. 4) GO TO 10 IF (NN .EQ. 8) GO TO 15 CALL TRICON(X,Y,S3,D,A3) GO TO 20 10 CALL QUACON(X,Y,S4,D,A4,H,TW,BC) GO TO 20 15 CALL LQCON(X,Y,S8,D,A8,H,TW,BC) 20 CONTINUE RETURN END Q *************************************************************** SUBROUTINE QUACON(X,Y,S,D,A4,H,TW,BC) C C CALCULATES CONDUCTANCE MATRIX FOR RECTANGULAR ELEMENT C IMPLICIT REAL*8(A-H,0-Z) DIMENSION S(4,4),X(4),Y(4),D(2,2),B(2,4),BT(4,2),BTD(4,2) 1,AJ(2,2),XI(2),AI(2,2),ANU(4),ANV(4),ST(4,4) 2,SP(4,4),A4(4),AN(1,4),ANT(4,1) DATA XI/-0.577350269l89626D0,+0.577350269189626D0/ C C BEGIN DO-LOOPS FOR NUMERICAL INTEGRATION 104 CALL DGSET(S,4,4,4,0.D0) DO 26 1=1,2 DO 26 J=1,2 CALL DGSET(AJ,2,2,2,0.D0) U = XI(I) V = XI(J) CALCULATES DERIVATIVES OF SHAPE FUNCTIONS ANU(1) = -(1.D0-V)*0.25D0 ANU(2) = -ANU(1) ANU(3) = (1.D0+V)*0.25D0 ANU(4) = - ANU(3) ANV(1) = - (1.D0 - U) *0.25D0 ANV(2) = - (1.D0 +U) * 0.25D0 ANV(3) = -ANV(2) ANV(4) = -ANV(1) CALCULATE JACOBIAN DO 6 K=1,4 AJ(1 , 1 ) = AJ(1 , 1) A J ( 1 , 2) = A J (1 , 2) A J (2 , 1 ) = A J (2 , 1 ) AJ(2,2) = AJ(2,2) 6 CONTINUE + ANU(K)*X(K) + ANU(K) * Y(K) + ANV(K) * X(K) + ANV(K) * Y(K) - CALCULATE DETERMINANT AND INVERT AJ DET = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) AI (1,1) = AJ(2,2)/DET AI(1,2) = - AJ(1,2)/DET AI(2,1) = -AJ(2,1)/DET AI(2,2) = AJ(1,1)/DET CALCULATE B DO 8 K=1,4 B(1,K) = AI(1,1)*ANU(K) +AI(1,2)*ANV(K) 8 B(2,K) = AI(2,1)*ANU(K) + AI(2,2)*ANV(K) CALCULATE BT*D*B CALL DGTRAN(B,BT,2,4,2,4) CALL DGMULT(BT,D,BTD,4,2,2,4,2,4) CALL DGMULT(BTD,B,ST,4,2,4,4,2,4) MULTIPLY BY DET AND WEIGHT DO 10 K=1,4 DO 10 L=1,4 10 S(K,L) = S(K,L) + ST(K,L)*DET 26 CONTINUE 105 C C C ICODE= BC CALL DGSET(SP,4,4,4,0.DO) DO 30 K=1,4 30 A4(K) = 0.D0 DO 50 1=1,2 CALL DGSET(AJ,2,2,2,0.D0) GO TO (34,33,32,31,50), ICODE 31 U = -1.DO V = XI(I ) ALINE = DABS(DSQRT((X(4)~X(1))**2 +(Y(4)-Y(1))**2)) GO TO 35 32 V = 1.DO U= XI(I) ALINE = DABS(DSQRT((X(3)~X(4))**2 + (Y(3)-Y(4))**2)) GO TO 35 33 U = 1.DO V = XI (I ) ALINE = DABS(DSQRT((X(2)-X(3))**2 + (Y(2)-Y(3))**2)) GO TO 35 34 V = -1.DO U = XI(I ) ALINE = DABS(DSQRT((X(1)-X(2))**2 + (Y(1)-Y(2))**2)) = (1.D0-U)*(1,D0-V)*.25D0 = (1.D0+U)*(1.D0-V) *.25D0 = (1.D0+U)*(1.D0+V)*.25D0 = (1.D0-U)*(1.D0+V)*.25D0 35 AN(1,1) AN(1,2) AN(1,3) AN(1,4) CALCULATE NT*N CALL DGTRAN(AN,ANT,1,4,1,4) CALL DGMULT(ANT,AN,SP,4,1,4,4,1,4) DO 40 K=1,4 DO 40 L=1,4 40 S(K,L) = S(K,L) + SP(K,L)*H*ALINE*.5D0 DO 45 K=1,4 45 A4(K) =A4(K) + AN(1,K)*H*TW*ALINE*.5D0 50 CONTINUE RETURN END Q *********** * * ************************************************** SUBROUTINE TRICON(X,Y,S,D,A3) C C CALCULATES CONDUCTANCE MATRIX FOR TRIANGULAR ELEMENT C IMPLICIT REAL*8(A-H,0-Z) DIMENSION S(3,3),A3(3),X(3),Y(3),D(2,2),B(2,3),BT(3,2), 1BTD(3,2) A= X(2)*Y(3) - Y(2)*X(3) + X(1)*(Y(2)-Y(3)) + 1Y(1)*(X(3)-X(2)) AINV = 1.DO/A B(1,1) = (Y(2) -Y(3)) *AINV B(1 ,2) = (Y(3)-Y(1)) * AINV B( 1 ,3) = (Y(1)-Y(2)) * AINV B(2, 1) = (X(3)-X(2)) * AINV 106 B ( 2 , 2 ) = ( X ( 1 ) - X ( 3 ) ) * A I N V B ( 2 , 3 ) = ( X ( 2 ) - X ( 1 ) ) * A I N V C C C A L C U L A T E B T * D * B C C A L L D G T R A N ( B , B T , 2 , 3 , 2 , 3 ) C A L L D G M U L T ( B T , D , B T D , 3 , 2 , 2 , 3 , 2 , 3 ) C A L L D G M U L T ( B T D , B , S , 3 , 2 , 3 , 3 , 2 , 3 ) D O 1 0 I = 1 , 3 D O 1 0 J = 1 , 3 1 0 S ( I , J ) = S ( I , J ) * 0 . 5 D 0 * A D O 2 0 1 = 1 , 3 2 0 A 3 ( I ) = 0 . D 0 R E T U R N E N D Q ************************************* S U B R O U T I N E L A Y O U T ( X , Y , I C O , N D I M , I T Y P E , N E , N N O D E S , B C T , B C C , H , T W 1 , C O T C , I H F , C O T A , I M A T ) C C T H I S R O U T I N E R E A D S I N D A T A F O R F I N I T E E L E M E N T H E A T C O N D U C T I O N C P R O B L E M C I M P L I C I T R E A L * 8 ( A - H , 0 - Z ) I N T E G E R B C C D I M E N S I O N X ( N D I M ) , I H F ( N D I M ) , Y ( N D I M ) , I C O ( N D I M , 8 ) , I T Y P E ( N D I M ) 1 B C T ( N D I M ) , B C C ( N D I M ) , I M A T ( N D I M ) R E A D ( 5 , 2 8 ) C O T C , C O T A 2 8 F O R M A T ( F 5 . 3 , F 8 . 6 ) W R I T E ( 6 , 2 9 ) C O T C 2 9 F O R M A T ( 1 X , ' C O - E F F . O F T H E R M A L C O N D U C T I V I T Y ( A L ) = ' , F 5 . 3 , / ) W R I T E ( 6 , 6 0 ) C O T A 6 0 F O R M A T ( 1 X , ' C O - E F F . O F T H E R M A L C O N D U C T I V I T Y ( A I R ) = ' , F 8 . 6 , / ) R E A D ( 5 , 3 0 ) H , T W 3 0 F O R M A T ( F 7 . 3 , F 5 . 1 ) W R I T E ( 6 , 3 2 ) H , T W 3 2 F O R M A T ( 1 X , ' H E A T T R A N . C O - E F F . = ' , F 7 . 3 , 3 X , ' W A T E R T E M P . = ' , 1 F 5 . 1 , / ) R E A D ( 5 , 4 0 ) N E , N N O D E S 4 0 F O R M A T ( 2 1 5 ) W R I T E ( 6 , 4 1 ) N E , N N O D E S 4 1 F O R M A T ( 1 X , ' N O . O F E L E M E N T S = ' , I 5 , 2 X , ' N O . O F N O D E S = ' , 1 5 , / ) W R I T E ( 6 , 4 4 ) 4 4 F O R M A T ( 6 X , ' X C 0 - 0 R D ' , 2 X , ' Y C O - O R D ' , 2 X , ' B O U N D A R Y 1 C O N D I T I O N S ' , / ) D O 1 0 1 = 1 , N N O D E S C C I F B C T D O E S N O T E X I S T , I T S H O U L D B E M A D E N E G A T I V E C R E A D ( 5 , 4 3 ) X ( I ) , Y ( I ) , B C T ( I ) 4 3 F O R M A T ( 3 F 7 . 3 ) W R I T E ( 6 , 4 5 ) I , X ( I ) , Y ( I ) , B C T ( I ) 4 5 F O R M A T ( 1 X , I 3 , 2 X , F 8 . 3 , 2 X , F 1 0 . 3 , 5 X , F 8 . 3 ) 1 0 C O N T I N U E 107 C C C C C C C C c c c c ITYPE IS 4 IF RECTANGULAR ; 3 IF TRIANGULAR READ(5,51) (ITYPE(I),I=1,NE) 51 FORMAT(20I2) WRITE(6,52) 52 FORMAT(/,/,1X,1 ELEMENT ',3X,'BOUNDARY',3X,'MATERIAL' 1,' NODE ') WRITE(6,53) 53 FORMAT(1X,' NO. ',3X,' CODE ',3X,' CODE ', 13X,' NUMBERS ',/) DO 11 1=1,NE IT = ITYPE(I) BCC > 1,2,3,4,5,: 1,2,3,4, REFER TO SIDES SUFFERING CONVECTION 5 IF THE ELEMENT SUFFERS NO CONVECTION IHF > 0 IF THERE IS NO NEED TO CALCULATE HEAT FLUX 1 OTHERWISE I MAT •> 0 FOR AIR 1 FOR COPPER C C C C C C C C C C C C C C C C C C READ(5,49) BCC(I),IHF(I),IMAT(I),(ICO(I,J),J=1,IT) 49 FORMAT(1515) WRITE(6,47) I,BCC(I),IMAT(I),(ICO(I,J),J=1,IT) 47 FORMAT(1X,I6,6X,I6,5X,I6,3X,8(I5,1,')) 11 CONTINUE RETURN END ****************************************** SUBROUTINE SETUP(NODES,NNODEL,A,NDIMA,B,NDIMB,PMAST,PELEM) PROGRAM TO SETUP MASTER STIFFNESS MATRIX FROM INDIVIDUAL FINITE ELEMENT STIFFNESS MATRICES. ALSO SETS UP LOAD VECTOR. NODES NNODEL NVAR A NDIMA B NDIMB PMAST PELEM INTEGER VECTOR CONTAINING ELEMENT NODE NUMBERS IN ORDER. NUMBER OF NODES PER ELEMENT NUMBER OF VARIABLES PER NODE MASTER STIFFNESS MATRIX WHICH MUST BE ZEROED BEFORE FIRST ENTRY TO SETUP FIRST DIMENSION OF A ELEMENT STIFFNESS MATRIX FIRST DIMENSION OF B MASTER LOAD VECTOR WHICH MUST BE ZEROED BEFORE FIRST ENTRY. ELEMENT LOAD VECTOR IMPLICIT REAL*8(A-H,0-Z) DIMENSION A(NDIMA,NDIMA),B(NDIMB,NDIMB),NODES(20), 1PMAST(NDIMA),PELEM(20) DO 1 1=1,NNODEL MM = NODES(I) 108 PMAST (MM) = PMAST (MM) +PELEM (I ) DO 1 J=1,NNODEL NN = NODES(J) 1 A(MM,NN)=A(MM,NN)+B(I,J) RETURN END Q ************************************* SUBROUTINE BOUND(I,BCT,AK,B,NDIM,NN) C C THIS SUBROUTINE MODIFIES MATRIX FOR FIXED TEMPERATURE C BOUNDARY CONDITIONS C IMPLICIT REAL*8(A-H,0-Z) DIMENSION BCT(NDIM),AK(NDIM,NDIM),B(NDIM) PERM = AK(I,I) DO 10 J=1,NN B(J) = B(J) -AK(J,I)*BCT(I) 10 AK(I,J) = 0.D0 DO 11 J=1,NN 11 AK(J,I) = 0.D0 AK(I,I) = PERM B(I) = AK(I,1)*BCT(I) RETURN END Q *************************************************************** SUBROUTINE HFLUX(IEL,NN,X,Y,T,D,NO) C C THIS ROUTINE CALCULATES HEAT FLUX C IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(NN),NO(8),HF(4),Y(NN),T(NN),HF4(4),D(2,2) 1,HF8(8) IF (NN.EQ. 3) GO TO 5 IF (NN .EQ. 8) GO TO 10 CALL FLUX4(X, Y, T,HF 4,D,HFI) WRITE(6,51) 51 FORMAT(1X,' ELEMENT NO. ',5X,' NODES ',' 1 HEAT FLOW ) WRITE(6,52) (NO(I),I=1,NN) 52 FORMAT( 1X, 1 6X, 4 (16 , 2X)) WRITE(6,53) IEL,(HF4(I),1=1,NN),HFI 53 FORMAT(1X,I6,1 OX,4(F6.2,2X),3X,F7.3) GO TO 20 5 CALL FLUX3(X,Y,T,HF 3,D) WRITE(6,54) 54 FORMAT(1X,' ELEMENT NO. ',5X,T HRAT FLUX ') WRITE(6,55) IEL,HF3 55 FORMAT(1X,16,12X,F7.3) GO TO 20 10 CALL FLUX8(X,Y,T,HF8,D,HFI) WRITE(6,56) 56 FORMAT(1X,' ELEMENT NO. ',5X,' NODES ',32X,' 1 HEAT FLOW' ) WRITE(6,57) (NO(l),1=1,NN) 109 57 FORMAT(1X,16X,8(I6,2X)) WRITE(6,58) IEL,(HF8(l),1=1,8),HFI 58 FORMAT(1X,I 6,10X,8(F6.2,2X),3X,F7.3) 20 CONTINUE RETURN END Q ******************************************** SUBROUTINE FLUX3(X,Y,T,HF,D) IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(3),Y(3),B(2,3),BT(2),T(3),BTD(2) A = X(2)*Y(3) - Y(2)*X(3) + X( 1 ) * (Y( 2)-Y ( 3) ) +Y(D* 1(X(3)-X(2)) AINV= 1.DO/A B(1 , 1) = (Y(2) -Y(3))*AINV B(1,2) = (Y(3) - Y(1))*AINV B(1,3) = (Y(1) - Y(2))*AINV B(2,1) = (X(3) -X(2))*AINV B(2,2) = (X(1) - X(3))*AINV B(2,3) = (X(2) - X(1))*AINV CALL DGMATV(B,T,BT,2,3,2) CALL DGMATV(D,BT,BTD,2,2,2) HF = BTD(1) RETURN END C *************************************************************** SUBROUTINE FLUX4(X,Y,T,HF,D,HFI) C C CALCULATES HEAT FLUX IN RECTANGULAR ELEMENT C IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(4),Y(4),T(4),AJ(2,2),AI(2,2),ANU(4),ANV(4) 1,XY(8),YX(8),YI(2),IC(4),B(2,4),XI(2),HF(4),BTD(2),BT(2), 2D(2 2) DATA YI/-.577350269189626D0, + .577350269189626D0/ DATA XI/-1.D0,+1.DO/ DATA IC/1,2,4,3/ IP = 0 DO 200 1=1,2 DO 200 J=1,2 CALL DGSET(AJ,2,2,2,0.D0) U = XI(I ) V = XI(J) ANU(1) = -(1.D0-V)*0.25D0 ANU(2) = -ANU(1) ANU(3) = (1.D0+V)*0.25D0 ANU(4) = - ANU(3) ANV(1) = - (1.D0 - U) *0.25D0 ANV(2) = - (1.D0 +U) * 0.25D0 ANV(3) = -ANV(2) ANV(4) = -ANV(1) C C CALCULATE JACOBIAN C DO 6 K=1,4 110 IA = IC(K) XY(K) = X(IA) YX(K) = Y(IA) AJ(1,1) = AJ(1,1) + ANU(K)*X(K) AJ(1,2) = AJ(1,2) + ANU(K) * Y(K) AJ(2,1) = AJ(2,1) + ANV(K) * X(K) AJ(2,2) = AJ(2,2) + ANV(K) * Y(K) 6 CONTINUE C C CALCULATE DETERMINANT AND INVERT AJ C DET = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) AI ( 1 , 1 ) = AJ(2,2)/DET AI (1 ,2) = - AJ(1,2)/DET AI (2, 1 ) = -AJ(2,1)/DET AI (2,2) = AJ(1,1)/DET DO 100 K=1,4 B(1,K) = AI(1,1)*ANU(K) + AI(1,2)*ANV(K) 100 B(2,K) = AI(2,1)*ANU(K) + AI(2,2)*ANV(K) CALL DGMATV(B,T,BT,2,4,2) CALL DGMATV(D,BT,BTD,2,2,2) IP = IP + 1 IK = IC(IP) HF(IK) = BTD(1) 200 CONTINUE HFI = 0.D0 ALINE2 = DABS(DSQRT((X(2) - X(3))**2 + (Y(2) - Y(3))**2)) DO 350 1=1,2 CALL DGSET(AJ,2,2,2,0.D0) U = 1.DO V = YI(I ) ANU(1) = -(1.D0-V)*0.25D0 ANU(2) = -ANU(1) ANU(3) = (1.D0+V)*0.25D0 ANU(4) = " ANU(3) ANV(1) = - (1.D0 - U) *0.25D0 ANV(2) = - (1.D0 +U) * 0.25D0 ANV(3) = -ANV(2) ANV(4) = -ANV(1) C C CALCULATE JACOBIAN C DO 250 K=1,4 AJ(1,1) = AJ(1,1) + ANU(K)*X(K) AJ(1,2) = AJ(1,2) + ANU(K) * Y(K) AJ(2,1) = AJ(2,1) + ANV(K) * X(K) AJ(2,2) = AJ(2,2) + ANV(K) * Y(K) 250 CONTINUE C C CALCULATE DETERMINANT AND INVERT AJ C DET = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) AI ( 1 , 1 ) = AJ(2,2)/DET AI (1 ,2) = - AJ(1,2)/DET 111 A l(2,1) = -AJ(2,1)/DET Al(2 , 2) = AJ(1 ,1)/DET DO 300 K=1,4 B(1,K) = Al(1,1)*ANU(K) + Al(1,2)*ANV(K) 300 B(2,K) = Al(2,1)*ANU(K) + Al(2,2)*ANV(K) CALL DGMATV(B,T,BT,2,4,2) CALL DGMATV(D,BT,BTD,2,2,2) HFI = HFI + BTD(1)*ALINE2*.5D0 350 CONTINUE RETURN END C ******************************************** SUBROUTINE LQCON(X,Y,S,D,A8,H,TW,BC) C C CALCULATES CONDUCTANCE MATRIX FOR 8 NODE C ISOPARAMETRIC ELEMENT. C IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(8),Y(8),S(8,8),D(2,2),A8(8),W(3),XI(3), 1AJ(2,2),AI(2,2),ANU(8),ANV(8),AN(1,8),ANT(8,1),ST(8,8), 2SP(8,8),BT(8,2),BTD(8,2),B(2,8) DATA W/0.5555555555556D0,0.8888888888889D0,0.5555555555556 1D0/ DATA XI/-0.774596669241483D0,0.DO,0.774596669241483D0/ C C ZERO ARRAYS C CALL DGSET(S,8,8,8,0.D0) C C BEGIN DO-LOOPS FOR NUMERICAL INTEGRATION C DO 26 1=1,3 DO 26 J=1,3 CALL DGSET(AJ,2,2,2,0.D0) U=XI(I) V= XI(J) C C CALCULATE DERIVATIVES FOR SHAPE FUNCTIONS C ANU(1) = (1.D0-V)*(2.D0*U+V)*.25D0 ANU(2) = (1.D0-V)*(2.D0*U-V)*.25D0 ANU(3) = (1.D0+V)*(2.D0*U+V)*.25D0 ANU(4) = (1.D0+V)*(2.D0*U-V)*.25D0 ANU(5) = -U*(1.D0-V) ANU(6) = (1.D0-V*V)*.5D0 ANU(7) = -U*(1.D0+V) ANU(8) = -(1.D0-V*V)*.5D0 ANV(1) = (1.D0-U)*(U+2.D0*V)*.25D0 ANV(2) = (1.D0+U)*(2.D0*V-U)*.25D0 • ANV(3) = (1.D0+U)*(2.D0*V+U)*.25D0 ANV(4) = (1.D0-U)*(2.D0*V-U)*.25D0 ANV(5) = -(1.D0-U*U)*.5D0 ANV(6) = -V*(1.D0+U) ANV(7) = (1.D0-U*U)*.5D0 1 12 - ANV(8) = -V*(1.DO-U) C C CALCULATE JACOBIAN C DO 6 K=1,8 AJ(1,1) =AJ(1,1) +ANU(K)*X(K) AJ(1,2) = AJ(1,2) + ANU(K)*Y(K) AJ(2,1) = AJ(2,1) + ANV(K) * X(K) 6 AJ(2,2) = AJ(2,2) + ANV(K)*Y(K) C C CALCULATE DET AND INVERT AJ C DET = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) Al ( 1 , 1 ) = AJ(2,2)/DET Al(1 ,2) = -AJ(1,2)/DET Al(2,1) = -AJ(2,1)/DET AI(2,2) =AJ(1,1)/DET C C CASLCULATE CONDUCTANCE MATRIX C DO 8 K=1,8 B(1,K) = Al(1,1)* ANU(K) +AI(1,2)* ANV(K) 8 B(2,K) = AI(2,1) *ANU(K) + AI(2,2) * ANV(K) C C MULTIPLY BT*D*B C CALL DGTRAN(B,BT,2,8,2,8) CALL DGMULT(BT,D,BTD,8,2,2,8,2,8) CALL DGMULT(BTD,B,ST,8,2,8,8,2,8) C C MULTIPLY RESULT BY WEIGTHS AND DET. OF J AND C ADD INTO CONDUCTANCE MATRIX C DO 11 K=1,8 DO 11 L=1,8 11 S(K,L) = S(K,L) +W(I)*W(J)*ST(K,L)*DET 26 CONTINUE ICODE = BC CALL DGSET(SP,8,8,8,0.D0) DO 30 K=1,8 30 A8(K) = 0.D0 DO 50 1=1,3 CALL DGSET(AJ,2,2,2,0.D0) GO TO (34,33,32,31,50), ICODE 31 U= -1.DO V= XI(I) GO TO 35 32 V= 1.DO U= XI(I) GO TO 35 33 U = 1.D0 V= XI(I) GO TO 35 34 V= -1.DO 1 13 U= XI(I ) 35 AN(1,1) = -(1.DO-U)*(1.DO-V)*(1.DO+U+V)*.25D0 AN(1,2) = -(1.DO+U)*(1.DO-V)*(1.DO-U+V)*.25D0 AN(1,3) = -(1.DO+U)*(1.DO+V)*(1.D0-U-V)*.25D0 AN(1,4) = -(1.DO-U)*(1.DO+V)*(1.DO+U-V)*.25D0 AN(1,5) = (1.D0-U**2)*(1.D0-V)*.5D0 AN(1,6) = (1.D0-V**2)*(1.D0+U)*.5D0 AN(1,7) = (1.DO-U**2)*(1.D0+V)*.5D0 AN(1,8) = (1.D0-V**2)*(1,D0-U)*.5D0 C C CALCULATE DERIVATIVES FOR SHAPE FUNCTIONS C ANU(1) = (1.D0-V)*(2.D0*U+V)*.25D0 ANU(2) = (1.D0-V)*(2.D0*U-V)*.25D0 ANU(3) = (1.D0+V)*(2.D0*U+V)*.25D0 ANU(4) = (1.D0+V)*(2.D0*U-V)*.25D0 ANU(5) = -U*(1.DO-V) ANU(6) = (1.D0-V*V)*.5D0 ANU(7) = -U*(1.DO+V) ANU(8) = -(1.D0-V*V)*.5D0 ANV(1) = (1.D0-U)*(U+2.D0*V)*.25D0 ANV(2) = (1.D0+U)*(2.D0*V-U)*.25D0 ANV(3) = (1.D0+U)*(2.D0*V+U)*.25D0 ANV(4) = (1.D0-U)*(2.D0*V-U)*.25D0 ANV(5) = -(1.D0-U*U)*.5D0 ANV(6) = -V*(1.DO+U) ANV(7) = (1.D0-U*U)*.5D0 ANV(8) = -V*(1.DO-U) C C CALCULATE JACOBIAN C DO 37 K=1,8 AJ(1,1) =AJ(1,1) +ANU(K)*X(K) AJ(1,2) = AJ(1,2) + ANU(K)*Y(K) AJ(2,1) = AJ(2,1) + ANV(K) * X(K) 37 AJ(2,2) = AJ(2,2) + ANV(K)*Y(K) IF (ICODE .EQ. 1) GO TO 38 IF (ICODE .EQ. 3) GO TO 38 ALINE = Aj(2,1) + AJ(2,2) GO TO 39 38 ALINE = AJ(1,1) +AJ(1,2) 39 CONTINUE C C CALCULATE NT*N C CALL DGTRAN(AN,ANT,1,8,1,8) CALL DGMULT(ANT,AN,SP,8,1,8,8,1,8) DO 40 K=1,8 DO 40 L=1,8 40 S(K,L) = S(K,L) +SP(K,L)*H*W(I)*ALINE DO 45 K=1,8 45 A8(K) = A8(K) + AN(1,K)*H*TW*W(I)*ALINE 50 CONTINUE RETURN 1 14 END C***************************************************** SUBROUTINE FLUX8(X,Y,T,HF,D,HFI) IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(8),Y(8),T(8),AJ(2,2),AI(2,2),ANU(8),ANV(8) 1, XY(8),YX(8),YI(3),IC(8),B(2,8),XI(3),HF(8),BTD(2),BT(2) 2, D(2,2),W(3) DATA YI/-0.774596669241483D0,0.DO,0.774596669241483D0/ DATA W/0.5555555555556DO,0.8888888888889DO,0.5555555555556 1D0/ DATA XI/-1.D0,0.D0,+1.DO/ DATA IC/1,5,3,8,6,4,7,3/ IP=0 DO 200 1=1,3 DO 200 J=1,3 CALL DGSET(AJ,2,2,2,0.D0) U= XI(I) V = XI(J) IF (U .EQ. 0) GO TO 2 GO TO 3 2 IF (V .EQ. 0) GO TO 200 3 CONTINUE C C CALCULATE DERIVATIVES FOR SHAPE FUNCTIONS C ANU(1) = (1.D0-V)*(2.D0*U+V)*.25D0 ANU(2) = (1.D0-V)*(2.D0*U-V)*.25D0 ANU(3) = (1.D0+V)*(2.D0*U+V)*.25D0 ANU(4) = (1.D0+V)*(2.D0*U-V)*.25D0 ANU(5) = -U*(1.D0-V) ANU(6) = (1.D0-V*V)*.5D0 ANU(7) = -U*(1.D0+V) ANU(8) = -(1.D0-V*V)*.5D0 ANV(1) = (1.D0-U)*(U+2.D0*V)*.25D0 ANV(2) = (1.D0+U)*(2.D0*V-U)*.25D0 ANV(3) = (1.D0+U)*(2.D0*V+U)*.25D0 ANV(4) = (1.D0-U)*(2.D0*V-U)*.25D0 ANV(5) = -(1.D0-U*U)*.5D0 ANV(6) = -V*(1.D0+U) ANV(7) = (1.D0-U*U)*.5D0 ANV(8) = -V*(1.D0-U) C C CALCULATE JACOBIAN C DO 6 K=1,8 IA = IC(K) XY(K) = X(IA) YX(K) = Y(IA) AJ(1,1) =AJ(1,1) +ANU(K)*X(K) AJ(1,2) = AJ(1,2) + ANU(K)*Y(K) AJ(2,1) = AJ(2,1) + ANV(K) * X(K) 6 AJ(2,2) = AJ(2,2) + ANV(K)*Y(K) C C CALCULATE DET AND INVERT AJ 115 C DET = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) Al(1,1) = AJ(2,2)/DET Al(1,2) = -AJ(1,2)/DET Al(2,1) = -AJ(2,1)/DET AI(2,2) =AJ(1,1)/DET DO 8 K=1,8 B(1,K) = AI(1,1)* ANU(K) +AI(1,2)* ANV(K) 8 B(2,K) = AI(2,1) *ANU(K) + AI(2,2) * ANV(K) CALL DGMATV(B,T,BT,2,8,2) CALL DGMATV(D,BT,BTD,2,2,2) IP = IP + 1 IK = IC(IP) HF(IK) = BTD(1) 200 CONTINUE HFI = 0.D0 DO 350 1=1,3 CALL DGSET(AJ,2,2,2,0.D0) U = 1.DO V= YI(I ) C C CALCULATE DERIVATIVES FOR SHAPE FUNCTIONS C ANU(1) = = (1.D0-V)*(2 .D0*U+V)* .25D0 ANU(2) = = (1.D0-V)*(2 .D0*U-V)* .25D0 ANU(3) = = (1.D0+V)*(2 .D0*U+V)* .25D0 ANU(4) = = (1.D0+V)*(2 .D0*U-V)* . 25D0 ANU(5) = = -U*(1.DO-V) ANU(6) = = (1.D0-V*V)* .5D0 ANU(7) = = -U*(1.DO+V) ANU(8) = = -(1.D0-V*V) *.5D0 ANV(1) = = (1.D0-U)*(U+2.D0*V)* .25D0 ANV(2) = = (1.D0+U)*(2 .D0*V-U)* .25D0 ANV(3) = = (1.D0+U)*(2 .D0*V+U)* .25D0 ANV(4) = = (1.D0-U)*(2 .D0*V-U)* .25D0 ANV(5) = = -(1.D0-U*U) *.5D0 ANV(6) = = -V*(1.DO+U) ANV(7) = = (1.D0-U*U)* .5D0 ANV(8) = = -V*(1.DO-U) c C CALCULATE JACOBIAN C DO 250 K=1,8 AJ(1,1) =AJ(1,1) +ANU(K)*X(K) AJ(1,2) = AJ(1,2) + ANU(K)*Y(K) AJ(2,1) = AJ(2,1) + ANV(K) * X(K) 250 AJ(2,2) = AJ(2,2) + ANV(K)*Y(K) C C CALCULATE DET AND INVERT AJ C DET = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) Al ( 1 , 1 ) = AJ(2,2)/DET Al(1,2) = -AJ(1,2)/DET Al(2,1) = -AJ(2,1)/DET 116 AI(2,2) =AJ(1,1)/DET DO 300 K=1,8 B(1,K) = AI(1,1)* ANU(K) +AI(1,2)* ANV(K) 300 B(2,K) = AI(2,1) *ANU(K) + Al(2,2) * ANV(K) CALL DGMATV(B,T,BT,2,8,2) CALL DGMATV(D,BT,BTD,2,2,2) HFI = HFI + BTD(1)*W(I)*(AJ(1,1)+AJ(1,2)) 350 CONTINUE RETURN END 1 17 APPENDIX E ~ F.E. FORMULATION FOR STEADY HEAT TRANSFER IN A PLANE a. The Fourier heat conduction equation for steady heat transfer in a plane may be written as [ Q x q y l T = " E k ] [dT/dx 9T/9y] T where [k] = k x k x y k k For homogeneous isotropic materials k x y = 0 and k x Then k(9 2T/9x 2) + k(9 2T/9y 2) = 0. A functional for th i s problem i s = k, n = 0.5 ff({9T/9x 9T/9y} [k] {9T/9x 9T/9y}T dA + fq* T dS + J(h/2)(T 2 - 2TT»)dS where q* - prescribed heat flux on boundary SN . h - heat transfer c o e f f i c i e n t on boundary Sc . and T» - the attendant bulk temperature of the surrounding medium. On minimisation, the matrix equation [K ] {T} = {P^} + {pr } i s obtained where {T} i s the solution for the temperature d i s t r i b u t i o n and [K ] i s the conductance matrix given by /[B ] T [k] [B] dA ( The s p a t i a l f i e l d for temperature in an element i s [N ] {T } where [N] i s the shape function matrix. Then {9T/9x 9T/9y} = [B ] {T } where [B ] = {9 /9x 9 /9y} [N ]) {P^} = /[N ] T q* dS 118 1 t =-1 N = I N, N 2 N 3 N 4] where N, = 1/4 ( 1 - s ) ( 1 - t ) N2= 1/4 (1 * s » 1-t) N 3 = 1/4 (1*s)(1*t) N t = 1/4 (1- s)(1*t| Figure 51 - Shape functions-Rect. isoparametric element N = 2 T o 7 i o - / N 1 N 2 N 3 J where 1 X 2 Y 3 - Y 2 X 3 Y 2 -Y 3 X 3 - X 2 X 3 Y r Y 3 X 1 V Y 1 X 1 - X 3 L X 1 Y 2" Y 1 X 2 V Y 2 V X 1 Figure 52 - Shape functions-Triangular element 119 (P r } = J[N ] T h T«, dS the l a s t two equations representing the boundary conditions, b. The shape functions for a rectangular isoparametric element and those for a triangular element are given in f i g . 51 and 52. Ref. R.D.Cook, Concepts and Applications of F i n i t e Element Analysis , J.Wiley and Sons, 1981.