Electron Beam Melting and Refining of Metals: Computational Modeling and Optimization

Computational modeling offers an opportunity for a better understanding and investigation of thermal transfer mechanisms. It can be used for the optimization of the electron beam melting process and for obtaining new materials with improved characteristics that have many applications in the power industry, medicine, instrument engineering, electronics, etc. A time-dependent 3D axis-symmetrical heat model for simulation of thermal transfer in metal ingots solidified in a water-cooled crucible at electron beam melting and refining (EBMR) is developed. The model predicts the change in the temperature field in the casting ingot during the interaction of the beam with the material. A modified Pismen-Rekford numerical scheme to discretize the analytical model is developed. These equation systems, describing the thermal processes and main characteristics of the developed numerical method, are presented. In order to optimize the technological regimes, different criteria for better refinement and obtaining dendrite crystal structures are proposed. Analytical problems of mathematical optimization are formulated, discretized and heuristically solved by cluster methods. Using important for the practice simulation results, suggestions can be made for EBMR technology optimization. The proposed tool is important and useful for studying, control, optimization of EBMR process parameters and improving of the quality of the newly produced materials.


Introduction
Electron beam melting and refining (EBMR) in a vacuum using an intense electron beam is a widely used, ecologically-friendly method in special electro metallurgy for new materials fabrication: the production of pure metals and special alloys [1][2][3][4][5][6]. The obtained materials have improved chemical compositions, structures, properties, and have many applications in medicine, aerospace engineering, nuclear industry, instrument engineering, etc. The electron beam melting and refining of metals is accomplished in a vacuum chamber [the working vacuum pressure is about (5-8) × 10 −3 Pa] using electron beams as a heating source. The raw material is melted, refined and re-solidified in a water-cooled copper crucible, which is used due to the increasing demand for purity of the transformed materials ( Figure 1). The electrons fall on the front side of the feeding rod (raw material) and heat it. Then drops of molten metal fall in the crucible. The top surface of the molten metal in the water-cooled crucible (G1) is also heated by the electrons (Figure 1). By using a water-cooled pulling mechanism (Figure 1), the operator can withdraw the bottom of the growing formed pure metal ingot. This allows the liquid pool surface to be maintained at a constant level suitable to be viewed by the operator. This method provides a high refining level (high level of impurities removal-gases, metals and non-metals), chemical composition homogeneity and optimal structure of the cast pure ingots [1][2][3][4][5][6]. Despite development of the technology, there are many unsolved problems concerning mechanisms and relationships of macroscopic heat and mass transfer during heating free liquid surface with intensive energy flow, that still exist [3][4][5][6][7]. For the EBMR process and equipment optimization and the price of the produced pure materials, the questions concerning the improved energy efficiency are crucial. The answer of these questions depends on the detailed study of the heat transfer processes taking place in the zone of interaction beam-material, on the processes and factors limiting the geometry of the molten pool and the precise evaluation of the temperature field dynamics and behavior of the metals and their compounds during the EBMR. Due to difficulties in acquiring real time information about the processes in the molten pool [1,2], the successful application and optimization of EBMR depends also on the adequate mathematical modeling of the processes that allow for studying the influence of many technological parameters and the various limiting factors. Stationary and quasi-stationary heat models in this field are already published [8][9][10][11][12][13][14][15] and reveal the importance of some process characteristics. In [16] a finite element based mathematical model to predict the evolution of the temperature and stress during EBMR of solar-grade silicon ingots is also presented. Applying a heuristic approach the flatness of the liquid/solid boundary is examined and studied [17,18].
A time-dependent thermal mathematical model of the electron beam drip melting and refining process is developed and described in this paper. The model is a continuation and an extension of our quasi-steady-state heat model [10,12,15]. The boundary conditions on the ingot interfaces are formulated assuming different mechanisms of the heat transfer. The temperature variations of the thermo-physical parameters (heat capacity, thermal conductivity, etc.) are taken into account in the heat model. An optimization scheme, based on the heat model, is proposed and applied for quality improvement of new materials obtained by EBMR.

Main Equation
Let H be the ingot's height, R-the ingot's radius (Figure 1), and F -the maximal heating time that is chosen. The investigated temperature distributions are along cylindrical ingots (samples) and an assumption for angle symmetry is applied. Then in: T (r, z, t) denotes the temperature at time t at the points with height z and polar distance r ( Figure 2). In the inner points of Ω, the temperature distribution T (x) is described by the heat Equation (2) in cylindrical coordinates: where ρ [kg/m 3 ] is the density of the metal; C p [W · s/kg · K] is the heat capacity; λ [W/m · K] is the thermal conductivity; a [m 2 /s] is the thermal diffusivity and a = λ/(ρC p ) = 1/k. The last term in Equation (2) V a ∂T ∂z shows the casting, i.e., the heat added by the poured molten metal (from the melting raw material- Figure 1) into the crucible. It is given by the thermal energy transfer from the material moving with velocity V , coincident to the z-axis (Figures 1 and 2). When V = 0, the model describes the solidification of the cast ingot in the water-cooled crucible-"drip melting process" (Figure 1), while the case V = 0 describes the "disks melting method" (no material is added to the ingot by pouring) [2,10]. The thermo-physical parameters λ and C p of the investigated metals are modeled as functions of the temperature (by linear regression method) using experimental data [19][20][21][22]. Figure 2. Cylindrical metal ingot, T (r, z, t)-the temperature at the moment t in the ingot's point with coordinates (r, z).

Equations for the Boundary Conditions
The boundaries G1, G2, G3, G4 ( Figure 1) are described as follows: (3) The boundary Equations that correspond to these interfaces are: The heat transfer through G1 boundary (Equation (7)), where the radiation losses predominate, is described by the Stefan-Boltzmann law [3,12]. In Equation (7) α is the metal's emissivity, σ = 5.6704 × 10 −8 [J/s · m 2 · K 4 ] is the Stefan-Boltzmann constant. P surf (r, t) is the e-beam power density function and usually Gaussian-like distribution is used [2], defined by the beam diameter's value. The vapor losses are calculated by C p · W v · T (Equation (7)). The weight loss velocity is evaluated using experimental data for the material losses obtained during EBMR process; T room is the room temperature (300 K).
Equation (8) describes Newton's type of heat transfer through G2-the interface molten metal/water-cooled crucible side wall. The gradient ∂T 1 /∂r is equal to λ 2 (T water − T (x))/∆wall; λ 2 is the conductivity of the crucible's material (copper). When the liquid metal pool does not reach the interface G2 and/or its volume is small this boundary condition can be neglected (i.e., Q = 0, Figure 1). ∆ wall is the width of the crucible side wall and T water is the mean water temperature in the water-cooling system, (T water = 300 K).
The thermal transfer through the interface ingot/ vacuum (G3) is described by the Stefan-Boltzmann law and is represented by Equations (9) and (10) describes the heat transfer through the ingot bottom (G4)-area with an ideal heat contact. The gradient ∂T 2 /∂z is equal to λ 2 (T water − T (x))/∆under, where ∆under is the width of the water-cooled puller ( Figure 1).
The initial condition gives information about the temperature distribution at the moment t = 0 s. Usually, the temperature along the metal ingot at the first moment is a constant and is equal to the mean room temperature T room . This condition is described by Equation (11):

Heat Streams
The heat streams through the boundaries at the moment t are calculated as follows: 2.3.1. Energy Losses (Radiation and Vapor Losses) through G1: 2.3.2. Streams through the Boundaries G2, G3 and G4: 2.3.3. The Incoming Heat is Determined by the Heating Beam Energy Distribution P surf and by the Heat Added by the Poured Molten Drops P add : The additional heat due to the added liquid drops is calculated by: where R rod is the radius of the feeding rod ( Figure 1); q [W · s/kg] is the specific melting heat of the metal; T melt is the melting temperature of the investigated metal; C * p is an average of the heat capacity. The heat transfer processes during EBMR become stationary when the difference between the input heat and total energy losses through the boundaries becomes close to 0.

Numerical Method
A modified Pismen-Rekford method [23] is developed and applied for solving the problem Equations (2) and (7)(8)(9)(10)(11) and for calculating the temperature fields in the cast ingot. A regular net of points W h 1 ,h 2 ,τ in the domain Ω concerning Equations (2) and (7-11) is made: T n i,j denotes the approximated temperature in the point (r i , z j , t n ) of the grid. Additionally, T n is the matrix of M × N points, corresponding to the temperature distribution at the moment t n = nτ . For the purposes of the algorithm intermediate planes T n+1/2 (which are not included in the net W h 1 ,h 2 ,τ ) are built. Approximations T n of T (r, z, τ n), continuously starting from n = 0 to n = P , are made. Geometrically, "jumps" are continuously made on rectangles Ω n = {(r, z, τ n)|0 ≤ r ≤ R, 0 ≤ z ≤ H}, parallel to Ω 0 . For each n, two jumps are associated and also two mathematical problems (A) and (B) have to be solved: Here Λ 1 and Λ 2 are difference operators, corresponding to 1 r ∂ ∂r r ∂T ∂r and V a ∂T ∂z , respectively:  i,j } i=0,N along Ω n (j and n are fixed, Figure 3). For the problem (B) such a system consists of equations concerning the points {T n+1 i,j } i=0,M (i and n are fixed, Figure 3). The obtained linear systems are three-diagonal and are solved via the Thomas method. The numerical scheme developed is absolutely stable and implicit in terms of τ due to the definition of problems (A) and (B) proposed in the non-stationary mathematical model. The thermal conductivity λ and the specific heat capacity C p for each investigated metal are implemented in the presented model as functions of the temperature (see 2.1, Table 1). When solving (A) and (B), temperature distribution from the previous iteration is used for evaluation of λ and C p .

Optimization Model and Different Criteria
The knowledge of the geometry of the crystallization front (liquid/solid boundary) is very important for studying and optimizing the quality of the obtained pure metal after the EBMR process. The flatness of the liquid/solid contour is directly connected to the quality of the structure of the obtained metal. The flat crystallization front permits the formation of vertical dendrite structure and uniform impurities' displacement toward the ingot top surface (Figures 4-6). When the liquid/solid boundary is deep in the center part of the metal block and shallow at the periphery of the cast ingot, a non-uniform structure along the ingot radius will be formed. The flatness of the temperature lines in a vertical cross-section of the cylindrical metal ingot for a fixed moment of time depends on the values of ∂T ∂r . For fixed moment of the heating t = t f and height z = z f , T (r; z f ; t f ) is a decreasing function of the variable r and ∂T ∂r ≤ 0. Our aim is to maximize the flatness of the contour liquid metal/ solid metal.

Different Criteria
In this paragraph we discuss different criteria aiming to obtain maximum flatness of the crystallization front shape. One approach is to minimize the average of − ∂T ∂r over an area that includes the liquid/ solid boundary. Then, the criterion is: where T is the temperature field determined by the solution of Equations (2) and (7-11) and µ(r, z, t) is a weight function. The control variables are the casting velocity V , the beam power P b and the beam radius r b · P b is connected with the heating beam energy distribution P surf . For a fixed heating time t f , the beam energy distribution P surf (r, t f ) is a Gaussian-like function and: Minimization of Equation (25) can be realized over an area D that includes the crystallization front (the molten pool contour, Figure 4a): In a more complicated case D (Figure 4b) can be defined as: Another approach is to minimize the mean curvature of the curve Γ(t f ), corresponding to the liquid/ solid boundary ( Figure 4) at a fixed moment of heating t f : where: or to minimize the mean curvature over a heating time interval: For Γ(t), s is a natural parameter for the curve and µ 1 , µ 2 are weight functions. The weight functions in Equations (25) and (31) can be defined in different ways taking into account the peculiarities of the electron beam melting process of the investigated metal. The used weight function in the paper is constant.

Discretization of the Criteria
The proposed criteria for achieving flatness of the crystallization front shape can be discretized synchronically to the proposed numerical scheme for the heat model. For the criterion (25) and D, described by Equation (27) or Equation (28), the multidimensional trapezoidal rule can be used. In the case of Equation (27) lets denote: where {T n i,j } n=1,P i=1,N ,j=1,M is the discrete temperature field calculated by the modified Pismen-Rekfort scheme. The subnet of W h 1 ,h 2 ,τ which is used to discretize the criterion is: In Equation (33), N i corresponds to R i , M i corresponds to H i , and P i corresponds to F i , i = 1, 2. Formula (34) presents the approximation of Equation (25): For discrete calculation of the criteria (29) and (31), the curve Γ is approximated by a set of points. Then the curve can be interpolated by linear regression method or some interpolation method.
After the criterion discretizing and choosing the optimization variables (the input power P b , the e-beam radius r b , the casting velocity V , etc.) for solving the defined optimization problem, an algorithm is needed. Heuristic methods (clustering optimization technique, genetic algorithms, etc.) are appropriate because the dependence of the criteria on the control variables is implicit.

Results and Discussion
The whole developed process tool of modeling, simulation and optimization of the process of EBMR of metals can be summarized as follows: 1. Theoretical mathematical modeling of the heating processes during EBMR of metals and alloys, describing the equations of the model; A corresponding computer program, based on the proposed time-dependent heat model, is developed for the study of the thermal processes at EBMR of metals. In this paper, results for EBMR of copper and hafnium are presented and discussed. The temperature variations of the thermal conductivity λ and the heat capacity C p for Cu and Hf are estimated using experimental data [19,21]. The obtained dependencies are presented in Table 1.    Table 2 at heating with different e-beam powers of 10 kW and 15 kW, r b = 12 mm and V = 0 mm/min. The obtained experimental results, investigated (after EMBR) by a metallographic etching method, are compared with simulation data, obtained by the presented heat model. The results about the geometry of the crystallization front shape-diameters d m and depths h m of the liquid metal pool are shown in Table 2. A good correspondence between calculated and experimentally obtained shapes of the crystallization front is observed. Using the developed simulation tool important data about the geometry of the liquid metal pool, heat streams through the boundaries, temperature fields in vertical and/or horizontal cross-sections of the cylindrical ingots during the EBMR process are obtained. Calculated and experimentally obtained crystallization front shapes are compared and a good correspondence is observed. Conclusions concerning the influence of a variety of regime parameters (such as e-beam radius, beam power, casting velocity, etc.) are made based on the simulation results [13,15,24]. They can be used for optimization of the heat streams according to the process requirements and for choosing proper process conditions. They are used in the developed optimization approach concerning the flatness of the crystallization front shape.
Using the developed tools and corresponding computer programs numerical experiments for EBMR of hafnium are made. The ingot's dimensions are H = 50 mm and R = 30 mm, the total heating time is 6 min, V = 0 mm/min. The radius of the e-beam is 18 mm and the beam power is 18, 20, 22, 24 kW. Simulation results for the heat transfer processes are investigated. It is observed that for each beam power the molten pool is mostly flat at the moment before the liquid pool reaches G2 boundary for a first time (Table 3). For all the beam powers the molten pool periodically reaches and withdraws G2. Table 3. First moment of contact between the molten pool and G2 boundary. Using the obtained simulation results, one dimensional cluster optimization technique for EBMR of Hf is performed. The control variable is the beam power P b . The chosen criterion is Equation (25) (Figures 4 and 5). The results obtained by the one dimensional cluster optimization show that for all three different areas D, the flattest molten pool (according to the criterion (25)) is observed for P b = 18 kW (Figure 6). So a possible technological suggestion is that for the investigated technological regimes a lower e-beam power (18 kW) has to be kept.
The cluster optimization technique is applied to optimize the technological regime for electron beam melting of titanium. The criterion (25) is minimized for P b = 11 kW [25]. This result is compared and coincides to the value obtained by experimental data and statistical method [26]. An important advantage of the optimization scheme proposed is that experimental data for chemical analysis of the impurities' concentrations for different technological regimes are not needed. The developed optimization approach can be applied to control and suggest proper technological parameters for electron beam melting of different metals for improving the quality of the obtained new materials.

Conclusions
A time-dependent thermal mathematical model, corresponding numerical method and computer program for simulation of heat transfer in metal ingots during EBMR are developed and presented. The heat model is 3D axis-symmetric and different heat transfer mechanisms through the boundaries are assumed. The developed numerical scheme is absolutely stable and implicit in terms of the time. The non-stationary heat model is programmed and corresponding computer software based on the model is also developed. The mathematical model allows acquiring important for the practice data and dependencies, which are otherwise difficult to obtain through EBMR experimental study (such as liquid pool geometry, energy losses, temperature distributions in metal ingots, etc.). The developed tool is employed to predict the evolution of the temperature, streams and crystallization front shape during electron beam melting and refining of different metals (hafnium, copper). Calculated results are compared to experimental data and good correspondence is observed. In order to optimize the EBMR process, different criteria, concerning the flatness of the crystallization front shape, are proposed and discussed. Analytical optimization problems are formulated, discretized using the heat model and solved by a cluster optimization method. The developed computational modeling and optimization tools and results obtained are important and give opportunity for better understanding, studying and optimizing the EBMR process and quality improvement of the purified materials produced by this expensive modern technology.