THE ADJOINT METHOD FOR THE DESIGN OF DIRECTIONAL BINARY ALLOY SOLIDIFICATION PROCESSES IN THE PRESENCE OF A STRONG MAGNETIC FIELD A Dissertation Presented to the Faculty of the Graduate School of Cornell University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy by Rajiv Sampath January 2001 c Rajiv Sampath 2001 ALL RIGHTS RESERVED THE ADJOINT METHOD FOR THE DESIGN OF DIRECTIONAL BINARY ALLOY SOLIDIFICATION PROCESSES IN THE PRESENCE OF A STRONG MAGNETIC FIELD Rajiv Sampath, Ph.D. Cornell University 2001 An adjoint method formulation and numerical implementation is presented for a specific class of inverse design solidification problems. In particular, a method is developed to calculate the thermal boundary conditions on the mold walls for directional solidification processes under the influence of an external magnetic field such that a desired stable solid-liquid interface growth is achieved throughout solidification. Achieving a stable interface growth has profound implications on the obtained cast microstructures and is directly related to the quality of the solidified product. The considered inverse solidification design problem belongs to a typical class of inverse heat transfer problems, in which, incomplete conditions are provided on one part of the boundary, whereas over-specified conditions are given on another part of the boundary. The design problem is mathematically posed as a functional optimization problem. The unknowns of the design problem are thermal boundary conditions on part of the mold walls. The cost functional is defined so as to represent the deviation of the freezing interface thermal conditions from local thermodynamic equilibrium. A continuum adjoint based method for gradient computations coupled with a conjugate-gradient optimization solver is employed to solve the inverse problem. The design method is developed in two stages. First, the adjoint method is formulated for an inverse magneto-convection problem in a fixed domain with convection driven by buoyancy effects as well as a Lorentz force generated due to the applied magnetic field. The developed methods are demonstrated using various examples in which the exact solution to the inverse problem is known a priori. The examples demonstrate the accuracy and convergence behavior of the method in the framework of the conjugate gradient algorithm. The need for regularization is identified in one of the examples with uniformly distributed random errors in the input/measured temperature data where a H 1 regularized formulation is introduced in order to obtain stable solutions. The method is shown to be very robust and to work well for various problems including 3D applications. Secondly, the above developed adjoint technique is applied to the design of two typical solidification design problems. A directional binary alloy solidification process is examined in which melt convection is induced due to the combined action of buoyancy as well as Lorentz forces due to an external magnetic field. The goal of the inverse design problem is to identify the thermal conditions on the mold walls so that a desired stable flat-interface growth is realized throughout solidification. The above method is then extended to the design of a directional solidification process of a near-eutectic binary alloy in which convection is driven by the coupled action of buoyancy, thermocapillary and electromagnetic convection. Finally, the thesis concludes with a discussion of possible extensions to the proposed method. Biographical Sketch The author was born in Madras, India on January 27, 1975. After completing his high school education from D. A. V. Higher Secondary School, the author was admitted into the bachelor's program at the Indian Institute of Technology, Madras in 1992, from where he received his B. Tech degree in Mechanical Engineering in July, 1996. In August 1996, the author was admitted into the Ph.D. program at the Sibley School of Mechanical and Aerospace Engineering, Cornell University. In April 1999 the author passed the admission to candidacy examination and was awarded a special Masters degree in May 1999. iii To my parents Malini Sampath and V. S. Sampath iv Acknowledgements I would like to thank my thesis advisor, Professor Nicholas J. Zabaras, for his support and guidance over the last 5 years. I would like to thank Professors Kenneth E. Torrance, Thomas Coleman and Subrata Mukherjee for serving on my special committee. I would like to also thank the Sibley School of Mechanical and Aerospace Engineering for supporting me with a teaching assistantship for part of my stay at Cornell. Partial support from ALCOA Laboratories and a research grant from the NASA microgravity materials science program (NRA-98-HEDS-05) is also gratefully acknowledged. The computing support and facilities provided by the Cornell Theory Center are acknowledged. Finally, I would like to thank my family and friends who helped and supported me. v Table of Contents 1 Introduction 1.1 Background . . . . . . . . . . . 1.2 Related previous works . . . . . 1.3 Objectives of the present study 1.4 Organization of the thesis . . . 2 An 2.1 2.2 2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 4 7 10 11 11 13 20 20 20 29 29 31 37 40 44 44 50 54 56 60 63 64 67 67 70 76 85 86 Inverse Magneto-Convection Problem Introduction . . . . . . . . . . . . . . . . . . . . . . . . Definition of the inverse magneto-convection problem . A functional optimization framework . . . . . . . . . . 2.3.1 Governing equations for the sensitivity problem 2.3.2 Gradient calculations . . . . . . . . . . . . . . . 2.4 Numerical implementation . . . . . . . . . . . . . . . . 2.4.1 The conjugate gradient algorithm . . . . . . . . 2.4.2 Implementation of the direct problem . . . . . . 2.4.3 Implementation of the sensitivity problem . . . 2.4.4 Implementation of the adjoint problem . . . . . 2.5 Numerical examples . . . . . . . . . . . . . . . . . . . . 2.5.1 Example 1 . . . . . . . . . . . . . . . . . . . . . 2.5.2 Example 2 . . . . . . . . . . . . . . . . . . . . . 2.5.3 Example 3 . . . . . . . . . . . . . . . . . . . . . 2.5.4 Example 4 . . . . . . . . . . . . . . . . . . . . . 2.5.5 Example 5 . . . . . . . . . . . . . . . . . . . . . 2.5.6 Computational issues . . . . . . . . . . . . . . . 2.6 Summary and conclusions . . . . . . . . . . . . . . . . 3 Inverse Design of Directional Binary Alloy Solidification Presence of a Strong Magnetic Field 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Direct alloy solidification in the presence of a magnetic field . 3.2.1 Simulation of a reference solidification problem . . . . . 3.3 Definition of the inverse alloy solidification problem . . . . . . 3.3.1 Enforcement of the constitutional stability condition . vi in the . . . . . . . . . . . . . . . . . . . . 3.4 3.5 3.3.2 Adjoint method formulation . . . . . . . Numerical implementation and results . . . . . 3.4.1 Definition of the inverse problem . . . . 3.4.2 Validation of the inverse design solution Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 91 96 103 104 4 Inverse Design of Directional Eutectic Solidification in the Presence of Buoyancy, Thermocapillary and Electromagnetic Convection 109 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2 Direct binary alloy solidification problem with coupled thermocapillary, buoyancy and electromagnetic melt convection . . . . . . . . . . 112 4.2.1 Mathematical model . . . . . . . . . . . . . . . . . . . . . . . 114 4.2.2 Solution procedure . . . . . . . . . . . . . . . . . . . . . . . . 117 4.2.3 Reference binary eutectic alloy solidification problem . . . . . 119 4.3 Definition of the inverse alloy solidification problem . . . . . . . . . . 125 4.3.1 Enforcement of the constitutional stability condition . . . . . 126 4.3.2 Adjoint variable method . . . . . . . . . . . . . . . . . . . . . 128 4.4 Numerical design example . . . . . . . . . . . . . . . . . . . . . . . . 134 4.5 Validation of the inverse design solution . . . . . . . . . . . . . . . . 140 4.6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . 147 5 Suggestions for Future Studies A Time integration scheme for the direct flow equations 149 153 B Derivation of Adjoint Equations for Binary Alloy Solidification under the Influence of an External Magnetic Field 155 C Derivation of Adjoint Equations for Binary Eutectic Alloy Solidification in the Presence of Buoyancy, Thermocapillary and Electromagnetic Convection 158 vii List of Tables
2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 3.6 3.7 4.1 4.2 4.3 4.4 Direct problem to define (x, t; qo ), v(x, t; qo ) and (x, t; qo ) . . . . Sensitivity problem to define (x, t; qo , qo ), V (x, t; qo , qo ) and (x, t; qo , qo ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Adjoint problem to define (x, t; qo ), (x, t; qo ) and (x, t; qo ) . . . The conjugate gradient algorithm . . . . . . . . . . . . . . . . . . . Computational statistics for the various example problems . . . . . . 21 . . . . 22 30 32 65 Alloy solidification problem in the presence of a magnetic field. . . . 75 Thermo-physical properties of germanium and characteristic scales . 78 Dimensionless parameters governing the reference solidification system 80 Direct problem to define (x, t; qol ), v(x, t; qol ), c(x, t; qol ) and (x, t; qol ) 92 Sensitivity problem to define (x, t; qol , qol ), V (x, t; qol , qol ), C (x, t; qol , qol ) and (x, t; qol , qol ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Adjoint problem to define (x, t; qol ), (x, t; qol ), (x, t; qol ) and (x, t; qol ) 94 The conjugate gradient algorithm . . . . . . . . . . . . . . . . . . . . 95 Physical constants for the Sb - 8.6% Ge system . . . . . . . . . Dimensionless groups and their characteristic values . . . . . . Adjoint problem to define (x, t; qol ), (x, t; qol ) and (x, t; qol ) The nonlinear conjugate gradient algorithm . . . . . . . . . . . . . . . . . . . . . . . 120 121 133 135 viii List of Figures
1.1 Schematic summary of various solidification morphologies and their dependence on the interfacial heat flux G and velocity v (after, Kurz and Fisher [17]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Schematic of the inverse magneto-convection problem. . . . . . . . . Schematic of the system in Example 1. . . . . . . . . . . . . . . . . . Given heat flux history qI (y, t) on the left vertical wall of the cavity (Example 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 k Initial guess qo = 0 and intermediate fluxes qo (y, t) at iterations 1, 10, 20 (Example 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) The optimal heat flux distribution calculated at the 35th CGM k iteration (b) Variation of the cost functional J (qo ) and of the norm k of the gradient J (qo ) with CGM iteration counter k (Example 0 1, qo = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 k Initial guess qo = 1 - t4 and intermediate fluxes qo (y, t) at iterations 1, 5, 20 (Example 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) The optimal heat flux distribution calculated at the 60th CGM k iteration (b) Variation of the cost functional J (qo ) and of the norm k of the gradient J (qo ) with CGM iteration counter k (Example 0 4 1, qo = 1 - t ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . � Heat flux history qI (y, t) (on side x = 0) and qo (y, t) (on side x = 1) calculated by solving a direct magneto-convection problem (Example 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 Initial guessed heat flux qo (Example 2). . . . . . . . . . . . . . . . . k Intermediate heat fluxes qo (y, t) at iterations 1, 10, 20 and 40 (Example 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) The optimal heat flux distribution calculated at the 46th CGM k iteration (b) Variation of the cost functional J (qo ) and of the norm k of the gradient J (qo ) with CGM iteration counter k (Example 2). Finite element mesh and schematic of the system analyzed in Example 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Given heat flux history qI (y, t) on the left vertical wall of the domain (Example 3). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 15 45 46 47 2.1 2.2 2.3 2.4 2.5 47 49 2.6 2.7 50 2.8 2.9 2.10 2.11 51 52 53 53 54 55 2.12 2.13 ix 0 k 2.14 Initial heat flux qo (y, t) = 0.0 and intermediate heat fluxes qo (y, t) at iterations 1, 10 and 30 (Example 3). . . . . . . . . . . . . . . . . . 2.15 (a) The optimal heat flux distribution calculated at the 50th CGM k iteration (b) Variation of the cost functional J (qo ) and of the norm k of the gradient J (qo ) with CGM iteration counter k (Example 3). 2.16 Schematic of the three-dimensional inverse problem in Example 4. Note that all sides that are not shaded are insulated. . . . . . . . . . 0 2.17 Initial guess heat flux qo (0.5, z, t) = 1 - t and the intermediate fluxes k qo (0.5, z, t) at iterations 1, 5, 10 (Example 4). . . . . . . . . . . . . . 0 2.18 Initial guess heat flux qo (y, 0.5, t) = 1 - t and the intermediate fluxes k qo (y, 0.5, t) at iterations 1, 5, 10 (Example 4). . . . . . . . . . . . . . 2.19 The optimal heat flux distribution calculated at the 20th CGM iter� ation: (a) qo (y, 0.5, t), (b) qo (0.5, z, t) (Example 4). . . . . . . . . . . � k 2.20 The variation of the cost functional J (qo ) and of the norm of the k gradient ||J (qo )|| with CGM iteration counter k (Example 4). . . . 2.21 The calculated heat flux solution at the 30th iteration with fluctuating temperature data and no regularization (Example 5). . . . . . . . 2.22 The calculated heat flux solution at the 30th iteration with fluctuating temperature data and an H 1 type regularization (Example 5). . . 55 56 57 58 59 59 60 61 63 71 79 80 81 82 83 84 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 Schematic of the direct alloy solidification problem in the presence of a magnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stream function plots at dimensionless times = 0.9 and = 1.5 for the alloy solidification problem with no magnetic field. . . . . . . Stream function plots at dimensionless times = 0.9 and = 1.5 for the alloy solidification problem with an applied magnetic field. . . Variation of the flow intensity with increasing Hartmann number at dimensionless time = 0.5 predicted by a finite element calculation. Isopleths at dimensionless times = 0.9 and = 1.5 for the alloy solidification problem with no magnetic field. . . . . . . . . . . . . . Isopleths at dimensionless times = 0.9 and = 1.5 for the alloy solidification problem with an applied magnetic field. . . . . . . . . . The standard deviation s (t) for the case with and without an applied external magnetic field. . . . . . . . . . . . . . . . . . . . . . . Stability region plots for the cases with no magnetic field and with applied magnetic field in the direction of solidification ((y, t) < 0 represents the stable domain). . . . . . . . . . . . . . . . . . . . . . . 84 x 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.18 3.19 4.2 Schematic of the inverse problem to achieve a sharp interface growth with a desired front velocity v f during the solidification of a binary alloy in a cavity under the influence an external applied magnetic field Bo : (a) problem requirements and unknowns; (b) inverse magneto-thermo-solutal convection problem in the liquid domain. The parameter < 0 is used here to define the desired level of stability. Also note that G and Gc are the thermal and solutal gradients at the interface in the direction of the normal n. All equations shown here are in dimensional form. . . . . . . . . . . . . . . . . . . . . . The desired over stability function (y, t) and the reference front velocity vr (t) derived by modifying the corresponding quantities based on the diffusion driven (no-convection) reference direct problem. . . Objective function variation with CGM iterations. . . . . . . . . . k Flux solution qol (y, t) at iterations (a) 5, (b) 10, (c) 20 and (d) 40 (optimal). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Interface concentration c(y, t) from magneto-thermo-solutal convection problem on a prescribed variable domain: (a) history; (b) contours of isopleths in the (y, t) plane. . . . . . . . . . . . . . . . . . Optimal heat flux qos (y, t) at the solid mold wall os . . . . . . . . . c Contours of (y, t) n - m n corresponding to the optimal design problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison of the calculated interface position for the initial and optimal design problems at dimensionless time = 0.5 with the `target state'. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sample stream function plots at dimensionless times = 0.5 and = 1.5 for the optimal design problem. . . . . . . . . . . . . . . . Sample isotherms at dimensionless times = 0.5 and = 1.5 for the optimal design problem. . . . . . . . . . . . . . . . . . . . . . . . . Sample isopleths at dimensionless times = 0.5 and = 1.5 for the optimal design problem. . . . . . . . . . . . . . . . . . . . . . . . . Typical phase diagram for binary alloys, which indicates the phases ^ as functions of temperature T and the concentration of the solute, ^ . . . . . . . . . . . . . . . . . . . . . . . . which is denoted by C. Schematic of the binary alloy solidification problem in an open-boat configuration under the influence of an externally applied magnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Equilibrium phase diagram for the Sb-Ge system. . . . . . . . . . . Contours of the fluid stream function (a-c), isotherms (d-f) and isocomps (g-i) for the reference eutectic solidification simulation. Simulation results are shown at dimensionless times = 0.1, = 0.8 and = 1.5. Note that clockwise circulations are denoted by dashed lines while counter-clockwise recirculating cells are depicted by solid lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 . 98 . 99 . 101 . 102 . 103 . 104 . 105 . 105 . 106 . 106 . 112 4.3 4.4 4.6 . 114 . 119 . 124 xi 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 4.18 4.19 Contours plot of the solid composition for the reference solidification problem. The contours are labelled in wt% Germanium. . . . . . . . 125 Examination of the constitutional stability assumption at the solidliquid interface for the reference design problem. Dimensionless con^ tours of G/|v f | + m(CE - Co )/DL are displayed. Note that stable growth ( 0) is achieved only for a very short period during the entire solidification process. . . . . . . . . . . . . . . . . . . . 126 Schematic of the inverse problem to achieve a sharp interface moving with a desired front velocity v f during solidification of a binary alloy in a cavity with an applied magnetic field Bo : (a) problem requirements and unknowns; (b) inverse subproblem in the liquid domain. Note that stability is enforced using the a-priori determined criterion for the quantity G/|v f |, which is based on the constitutional stability criterion. All equations shown here are in dimensional form. 129 The desired interface velocity v f (t) and the heat flux G calculated through the solution of a "diffusion-based" problem. . . . . . . . . . 137 k Variation of the cost functional J (qo ) and of the norm of the gradient k ||J (qo )|| with conjugate gradient iterations. . . . . . . . . . . . . . . 138 k Flux distribution qol (y, t) obtained at intermediate conjugate gradient iterations k = 1, 2, 5, and 10, starting from an initial guess heat 0 flux qol (y, t) 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 Optimal flux distribution qol (y, t) on the liquid side mold boundary. . 139 � Optimal flux distribution qos (y, t) on the solid side mold boundary. . 141 � Standard deviation v (t) for the optimal and reference design problems.143 Comparison of the solid domain grids at intermediate times = 0.1 and = 0.4 demonstrating the prominent reduction in interface curvature due to the application of the optimally calculated heat fluxes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 Contours of the fluid stream function (a-c), isotherms (d-f) and composition (g-i) for the optimal design solidification simulation. Simulation results are shown at dimensionless times = 0.1, = 0.8 and = 1.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 Examination of the constitutional stability assumption at the solidliquid interface for the optimal design solidification problem. Dimen^ sionless contours of G/|v f | + m(CE - Co )/DL are displayed. Notice that stable growth is achieved for the entire duration of solidification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 Contour plot of the solid composition for the optimal solidification design problem. The contours are labelled in wt% Germanium. Note the vertical uniformity in the solid composition in comparison to the solute distribution patterns obtained for the reference eutectic solidification problem. . . . . . . . . . . . . . . . . . . . . . . . . . . 146 xii Chapter 1 Introduction
1.1 Background There have been extensive applications of inverse problem theory in the context of partial differential equations (PDEs) since the pioneering work of Hadamard [1]. In contrast to the well-studied and well-posed direct problems, inverse problems are "illposed" in nature. The ill-posedness in the system arises due to non-uniqueness/nonexistence of the solution or as an instability of the solution to small changes in the input data [2, 3]. Despite such mathematical difficulties, inverse problems have been under extensive investigation in the past few decades. Inverse problem theory has been applied to various branches of engineering, including heat transfer [4], signal/image processing [5], wave-propagation [6], oceanography [7], ground-water modeling [8], and many others. In the context of inverse heat transfer problems, the inverse heat conduction problem (IHCP) [9] has been the focus of attention due to its simplicity and widespread engineering applications. The majority of inverse problems examined in this area can be generally categorized as either inverse identification problems or inverse de- 1 2 sign/control problems. Typical identification problems of heat conduction include the reconstruction of the boundary heat flux conditions from temperature measurements within a conducting body, determining thermal properties (such as thermal conductivity or heat capacity of solids) as a function of temperature [10] and the calculation of distributed heat sources from the knowledge of the temperature field at some points within the domain and/or on the boundary [4]. The level of error in the measured input data is crucial to the accuracy of the results of the inverse analysis [11]. In contrast, the objectives of the inverse design heat conduction problem involves determining the necessary system characteristics (for example, the boundary temperature/heat flux conditions, location and intensity of heat sources within the domain, etc.) such that a desired state or effect is realized. Of the several inverse design heat conduction problems studied, the problem of calculating the temperature/heat flux conditions on the boundary of a domain is of significant interest. The solution to this class of inverse design problems involves determining the heat flux function belonging to an admissible function space which minimizes the value of an appropriate cost/objective functional. Various solution methods such as the time sequential method [9, 12], the whole time domain method [4], including several finite and infinite-dimensional optimization schemes were developed to calculate the unknown functions and thus solve the inverse design problem. A particularly popular approach is the infinite-dimensional iterative regularization scheme, also referred here as the adjoint method. This method allows the analytical calculation of the gradient of the objective function via the solution of an appropriately defined continuum adjoint problem. Adjoint methods have been applied successfully to several inverse heat transfer problems [4] as well as other 3 complex physical models [13]. One of the main advantages of the adjoint system method over other design optimization techniques, such as, the sensitivity equation method [14], is that the computational cost of evaluating the gradient of the objective function is approximately the same as the cost of the functional evaluation through the solution of the "direct problem". Due to this property, the adjoint equation approach is usually preferred in dealing with large scale inverse problems such as inverse thermal design problems. Adjoint methods have also been recently applied to the design of solidification processes. Solidification is an important phase change process with applications in materials processing, metallurgy, purification of metals, growth of pure crystals from melts and solutions, solidification of castings and ingots, food processing, etc. Several complicated physical mechanisms such as heat and mass transfer, thermodynamics, fluid flow, phase transformation kinetics, capillarity, etc., play an important role in the solidification process. Due to the combined action of these various driving mechanisms, solidification processes are characterized by the existence of a whole range of length and time scales. References [15], [16] and [17] provide an excellent overview of the physics and applications of solidification phenomena. Currently, exact solutions to the governing physical equations exist only in the form of similarity solutions such as the Neumann solution for Stefan problems [18, 19]. Several researchers have worked on developing mathematical and numerical modelling techniques for computer simulation of these processes and there has been tremendous progress in this area (see for example, [20]) and new methods are actively being investigated. In the present study inverse methods were developed for the optimal design of a class of directional solidification processes subjected to the effects of an externally 4 applied magnetic field. It is a standard practice in several solidification and crystal growth systems to introduce an external magnetic field to stabilize the flow and reduce macrosegregation. Here, we present a systematic and mathematically rigorous method to address the design of this class of directional solidification processes. 1.2 Related previous works The classical IHCP techniques have been extended by several researchers to solve inverse conduction-based phase change problems [21, 22, 23, 24, 25, 26]. These problems can be generally classified as inverse Stefan problems, where "Stefan problem" refers to a conduction-based moving boundary problem describing the phase-change of a pure substance. Several issues related to the mathematical theory of existence/uniqueness of the solution to the Stefan problem are addressed in [27]. Within the general category of inverse Stefan problems, three different types of inverse problems have been identified and solved. These include: � Calculation of the melting/solidification front location from temperature measurements within the solid domain [25, 28, 29]. � Reconstruction of boundary heat transfer conditions from temperature measurements within the solidifying domain [21]. � Inverse design solidification problem to control the boundary heating/cooling conditions such that a desired solid-liquid interface growth velocity is achieved [22, 23, 26, 30]. The third category of inverse problems is very closely related to the present study. Since the solid-liquid interface velocity/location is provided explicitly as a part of 5 the design objectives, in this class of inverse problems, the variable domain of the solid and liquid phases are explicitly known. Using this concept, Reference [30] addressed an inverse design problem for a directional solidification process with three boundaries insulated, and one boundary in which thermal boundary conditions need to be calculated such that the resulting growth velocity is as close as possible to the desired front velocity. This inverse solidification design problem can be separated into a direct problem in the melt domain and an inverse heat conduction problem in the solid phase. In the above work, the effects of buoyancy-driven convection was also included in the direct problem in the liquid domain. Through the solution of this coupled heat transfer and fluid flow problem, the liquid-side interface heat flux was evaluated. The solid-side interface flux is then calculated using the Stefan condition (interfacial thermal balance) and the desired freezing interface velocity. This calculated solid-side interfacial flux along with the thermodynamic equilibrium condition on the interface together constitute the over-specified input data required to solve the inverse heat conduction problem in the solid. A whole time domain adjoint method was employed in [30] to solve the above mentioned inverse heat conduction problem. In [31], Yang and Zabaras extended the work presented in [30] to address the design of directional solidification processes for pure materials such that a desired freezing front heat flux and growth velocity are achieved in the presence of melt convection. This problem was reduced to a inverse heat conduction problem in the solid and an inverse natural convection problem in the melt phase. A new methodology proposed by the same authors in [32] to solve inverse convection problems was used in combination with IHCP solution techniques to solve the above design problem. As in earlier works, a whole time domain adjoint method coupled with a 6 conjugate-gradient algorithm was used in solving the above inverse design problem. Reference [33] extended the above method to address the design of directional dilute binary alloy solidification processes. Solidification of a binary alloy differs in many respects from solidification of pure substances. Usually, the phase transition takes place over a range of temperatures, rather than at a discrete temperature as in pure substances. Furthermore, since the chemical components have different solubilities in the solid and liquid phases, during phase change a chemical species may be preferentially incorporated or rejected at the solidification front. Also, in many solidification situations, the solidification front is not always smooth and, under most practical conditions, a variety of microscopically complicated growth structures develop (e.g., cellular, dendritic or globulitic). This transition from planar to cellular/dendritic morphology is a consequence of the so called "morphological instability" [17] of the solidifying interface, leading to a significant depreciation in the quality of the final product. Reference [33] addressed the important problem of controlling the thermal boundary conditions for a specific class of directional binary alloy solidification processes, such that, upon application of the optimally designed fluxes, a stable planar solid-liquid interface growth is achieved. In all the above investigations, a moving/deforming finite element method has been employed for the solution of the field equations in the inverse design problem. This technique was originally developed for two-dimensional Stefan problems [34], and extended to phase change problems of pure substances with natural convection in the liquid domain [35]. As a "front tracking" method for moving boundary problems [36], the advantages of the deforming FEM include the high accuracy in front location determination and the ease in managing the deformed elements. Furthermore, for inverse design solidification problems with prescribed front veloc- 7 ity/location, the deforming FEM is most appropriate due to its flexibility in handling the decoupled inverse subproblems in the solid and liquid region as two independent moving boundary problems. 1.3 Objectives of the present study It is well understood that freezing front velocity (v) and heat flux (G) are very important solidification parameters [16, 17] that determine the crystallographic structure and scale of cast microstructures in the solidified product. Figure 1.1 shows the various microstructures which are obtained for a typical alloy, with a melting range of 50 K, when G and v are varied. Assuming an imposed unidirectional heat flow, the product Gv, equal to the cooling rate, determines the scale of microstructures obtained. On the other hand, the quantity G/v determines the type of growth morphology (e.g., dendritic, cellular). In solidification processing, the crystallographic morphology, scale of microstructures and grain orientation are very important parameters that are directly related to the strength and mechanical properties of the final solidified product. Hence, a proper control of these parameters is necessary to obtain an improved product. These factors influenced the choice of both the interface heat flux G and growth velocity v as part of the design objectives in the earlier studies. For the same reasons, these parameters play a similar role in the present investigation. Fluid flow in the liquid melt during solidification is another physical process that plays a very important role in nucleation, growth, coarsening phenomena and also affects micro- and macrosegregation in metallic alloys. Buoyancy-driven melt flow due to thermal/solutal gradients is a very common driving force for melt convection in almost all solidification processes. Melt flow due to surface tension gradients 8 100 V (mm/s) eq ua 1 cast xi ing ed ed de nd rit ic ie or 0.01 Directional Solidification de ed ce nt ie or llu la r nd 0.0001 0.1 10 V = G D / T o rit e nt . T = co ns t 1000 G (k/mm)
Figure 1.1: Schematic summary of various solidification morphologies and their dependence on the interfacial heat flux G and velocity v (after, Kurz and Fisher [17]) 9 (Marangoni convection) is also common in some solidification systems. External forcing agents such as rotation, electromagnetic fields, shear forces, etc., can also act as driving agents for melt flow in the liquid domain. In all the previous studies, only buoyancy-driven convection was considered. One of the main contributions of the present study is the extension of the developed methods to include the effects of coupled convective flow due to electromagnetic as well as surface-tension forces. In the first part of the current study, the methods addressed in [33] to design directional binary alloy solidification processes with a desired G/v state are extended to include the additional influence of an external magnetic field. Magnetic fields have been imposed in the design of several casting and solidification systems to damp melt convection and hence improve the cast microstructures. The particular design solidification problem addressed is an inverse problem in which the objective is to calculate the thermal boundary conditions such that a desired stable growth is achieved in the presence of coupled electromagnetic and buoyancy-driven convection in the melt. The inverse solidification problem is separable into an inverse heat conduction problem in the solid and an inverse magneto-convection problem in the liquid domain. A substantial amount of effort in this study is directed towards developing a systematic methodology for solving this inverse problem involving melt convection due to coupled effects of electromagnetic and buoyancy forces. The whole time domain adjoint method is considered appropriate for the solution of this class of problems. In the second part of this thesis, the developed methods are extended to the inverse design of two-dimensional binary eutectic solidification processes driven by the coupled action of buoyancy, thermocapillary and electromagnetic forces. This work represents a preliminary step towards inclusion of several coupled forcing agents driving fluid motion into the design methodology. 10 1.4 Organization of the thesis In the subsequent chapter (Chapter 2), the mathematical formulation and numerical algorithms for an inverse magneto-convection problem is very systematically and rigorously developed [37]. This methodology forms the basic structure for the following investigations. This inverse methodology is then tailored towards addressing the design of solidification processes. In particular, we address the following problems: � An inverse design method (Chapter 3) for the design of dilute binary alloy solidification processes in the presence of electromagnetic and buoyancy-driven convection in the melt such that a desired stable growth is achieved � An inverse design method (Chapter 4) to achieve desired growth conditions (i.e. desired interfacial flux G and growth velocity v) in a binary eutectic solidification process driven by the coupled action of buoyancy, thermocapillary and electromagnetic forces. Finally, the present work concludes with a list of suggestions for possible future studies. Chapter 2 An Inverse Magneto-Convection Problem
2.1 Introduction Considerable research interest has developed recently in solving inverse convection problems due to their importance and relevance to several practical engineering systems and processes. Cuvelier [38] pioneered such inverse transient convection problems from an optimal control viewpoint. More recently, Moutsoglou [39] investigated the problem of recovering the boundary heat flux distribution in a steady free convection problem. Gunzburger and colleagues [40, 41] have addressed some interesting design and optimal control inverse convection problems in recent years. They considered the problem of calculating the temperature in the inflow boundary in order to control temperature peaks along bounding surfaces of containers of fluid flows. In the definition of their problem, the flow field was assumed known and it was not dynamically coupled with the heat transfer analysis. Very recently, Zabaras and Yang [31, 32] proposed an inverse design methodology for the design of 11 12 free-convection systems (where the fluid flow analysis is coupled to the heat transfer calculations) and they also demonstrated the applications of these methods to the design of solidification processes for pure materials. In this chapter, the inverse free convection methodology proposed in [32] is extended to include the effects of an external magnetic field. In particular, a free convecting, incompressible, electrically conducting, viscous fluid in the presence of a strong, externally applied magnetic field is considered. The design problem of interest refers to the calculation/identification of the heat flux history in part of the boundary that results in a desired temperature distribution in another part of the boundary, where the heat flux is prescribed. The present problem can also be interpreted in the form of calculating the thermal conditions on one part of a liquid container that is not accessible to experimentation, using temperature and heat flux measurements on another part of the boundary. This problem has applications in the design of several fluid flow systems employing a strong magnetic field. Examples include the design of self-cooled liquid-metal blanket of a magnetically confined fusion reactor [42] and a number of applications in materials processing including electromagnetic stirring of continuous castings [43] and control of microstructures in semiconductor crystal growth [44]. The sequencing of the various sections in this chapter is as follows. The physics of the problem of interest is presented in Section 2.2 together with the governing equations. Also, the design problem is stated as an optimization problem in L2 with appropriate definitions of the cost functional and its gradient. Temperature, fluid velocity and electric potential sensitivity fields are also defined and the need for the calculation of the adjoint to the sensitivity fields is discussed. The inverse magneto-convection problem is posed as a functional optimization 13 problem in Section 2.3, where an expression for the exact L2 gradient is obtained based on the adjoint method. In Section 2.4, the conjugate gradient algorithm as is applied to the solution of the design problem is reviewed. Also, this section addresses the finite-element formulation and numerical solution of the direct, adjoint and sensitivity subproblems. A stabilized Galerkin formulation with SUPG (streamlineupwind/Petrov-Galerkin) and PSPG (pressure-stabilizing/Petrov-Galerkin) stabilizing terms is used for the solution of the flow equations. A SUPG formulation with weighting functions consistent with the fluid flow equations is used for the thermal problem. The solution of the electric potential equation is achieved through standard Galerkin techniques. In Section 2.5, the algorithm is applied to various example problems with known exact solutions. These include a three-dimensional inverse magneto-convection problem and an inverse magneto-convection problem on an arbitrary curved domain, which clearly demonstrate the applicability of the proposed methodology to the design/identification of complex fluid flow systems. In addition, an example problem with errors in the input data is also investigated. This example is used to demonstrate the need for regularization in problems involving large fluctuating measurement errors. The results obtained in all examples are very good; however, the method is shown to be computationally intensive. Finally, in Section 2.6, the current work is summarized. 2.2 Definition of the inverse magneto-convection problem Let be a closed bounded region in nsd , where nsd is the number of space dimen- sions, with a piecewise smooth boundary (Fig. 2.1). The region is occupied by an 14 incompressible electrically conducting fluid and is subject to an external magnetic field in a non-isothermal environment. The motion of this fluid is initially driven by temperature induced density gradients. Motion in the presence of a magnetic field will in turn give rise to a Lorentz or Laplace force, which acts on the fluid so that an extra body force term F appears in the Navier-Stokes equation. The Lorentz force term F in such a flow is given as follows: F = e E + J � B (2.1) where e is the electric charge density of the fluid, E = - the electric field intensity, the electric field potential, J the electric current density and B the magnetic field. On the other hand, the electric current density is governed by Ohm's law for a moving medium: J = e v + e (- + v � B) (2.2) where e is the electrical conductivity and v the fluid velocity vector. In addition to the applied magnetic field Bo , there is an induced magnetic field produced by the electric currents in the liquid metal. We assume in the following that the walls of the cavity are electric insulators and the magnetic Reynolds number Rem = P rm Red is sufficiently small that the induced magnetic field is negligible with respect to the imposed constant magnetic field Bo . Here, the magnetic Prandtl number is given by P rm = /m , where is the viscosity and m is the magnetic diffusivity. The dynamic Reynolds number is given by Red = U L/, where U and L are some characteristic velocity and length scale, respectively. As e is usually very small in liquid metals, the terms e E and e v in Eqs. (2.1) and (2.2), respectively, are neglected [45, 46]. In addition to Ohm's law, the electric current density J is governed by the 15 conservation of electric current. �J=0 (2.3) Equation (2.3) can be used to eliminate the current density J from the system of governing equations and to write all equations in terms of the electric field potential . The governing equation for this scalar potential is as follows: 2 = � (v � B) (2.4) Known g 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 n 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000n 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 1111111111111111 0000000000000000 n Unknown qo g h1 h0 I Bo Known q I Known m Known q1 Figure 2.1: Schematic of the inverse magneto-convection problem. Thus, under the influence of any external magnetic field, Eq. (2.4) gives the potential function (if the velocity is known) and this function along with the magnetic field determines the Lorentz force as given by Eq. (2.1). The remaining governing equations are introduced next to complete the definition of the problem. Let be the density, k the conductivity and ( = k/c) the thermal diffusivity of the liquid in . All the thermo-physical properties are assumed to be constant. The characteristic scale for time is taken as L2 / and for velocity as /L. The dimensionless temperature is taken as = (T - Tref )/(Tin - Tref ) where 16 T , Tin and Tref are the temperature, initial temperature and reference temperature, respectively. The characteristic scale for the electric potential is taken as Bo where Bo is the externally applied magnetic field. Since only dimensionless quantities will be used from here on, the symbol is also used to denote the dimensionless electric potential. The key dimensionless quantities are the Prandtl number (P r), the Rayleigh number (Ra) and the Hartmann number (Ha). They are defined as P r = /, Ra = g(Tin - Tref )L3 / and Ha =
e 1/2 Bo L , respectively, where is the thermal expansion coefficient and g is the gravity constant. The incompressible Navier-Stokes equations are used for modelling the flow induced by both density variations and Lorentz force. It is assumed that the flow is laminar and it has no viscous dissipation nor Joule heating effects. The conservation of momentum for the velocity field v(x, t) in (x, t) � [0, tmax ] is given by v + (v)v = � - P rRaeg + Ha2 P r [- + v � eB ] � eB , t where the constitutive equation defining the stress tensor is given as = -pI + 2P r(v), (v) = 1 [v + (v)T ]. 2 The incompressibility condition takes the form � v(x, t) = 0, (x, t) � [0, tmax ]. (2.8) (2.6) (2.7) (2.5) In the equations above, I is the second order unit tensor and eg and eB are the unit vectors in the direction of gravity and in the direction of the externally applied magnetic field, respectively. 17 The temperature field in �[0, tmax ] is governed by the following energy equation + v � = � . t (2.9) The electric potential field in � [0, tmax ] is governed by the following Poisson equation as described earlier 2 = � (v � eB ). The known initial conditions are the following: v(x, 0) = v i (x), and (x, 0) = i (x), x . (2.12) x , (2.11) (2.10) The initial potential function can be easily solved for when the velocity field is known. The fluid velocity, v, is assumed known on the boundary . The no-slip condition is used here, i.e. v(x, t) = 0, (x, t) � [0, tmax ]. (2.13) The electric potential function (x, t) is governed by insulating wall conditions on the boundary: (x, t) � n = 0, (x, t) � [0, tmax ]. (2.14) In the part h of the boundary , we assume that a heat flux boundary condition is applied, while in the remaining part of the boundary, g , a temperature boundary condition is considered, i.e. h g = and h g = (see Fig. 2.1). However, the distribution of the boundary heat flux on h0 h is not known (h0 h1 = h , h0 h1 = ). The known flux and temperature fields are given below: (x, t) = g (x, t), (x, t) g � [0, tmax ], (2.15) 18 q(x,