NUMERICAL SIMULATION OF SOLIDIFICATION AND MELTING PROBLEMS USING ANSYS FLUENT 16.2 BY SHUBHAM PAUL (10300712142) DEBORIT DE BISWAS(10300712109) WASIM SAJJAD(10300712153) RAKESH KUMAR JHA(10300612037)
Under the Supervision of Astt. Prof. Debraj Das
DEPARTMENT OF MECHANICAL ENGINEERING HALDIA INSTITUTE OF TECHNOLOGY HALDIA-721657 MAY, 2016
CERTIFICATE
This is to certify that the work contained in the thesis entitled Numerical Simulation of “
Solidification and Melting Problems using Ansys Fluent 16.2 by Shubham Paul ”
(University Roll no. 10300712142), Deborit De Biswas ( University Roll no . 10300712109), Wasim Sajjad ( University Roll no. 10300712153), Rakesh Kumar Jha (University Roll no. 10300612037) of the Department of Mechanical Engineering, Haldia
Institute of Technology in partial fulfillment of the the requirements for the award of Bachelor of Technology Degree in Mechanical Engineering during the academic session 2012-2016 is a bonafide record of thesis work carried out by them under my supervision and guidance. Neither of this report nor any part of it has been submitted for any an y degree or o r any academic acad emic award elsewhere.
…………………………………………… .. …………………………………………….. Counter signed by Head of the Department
………………………. ………………………. Mr. Debraj Das (Thesis Advisor)
Acknowledgement We would like to take this golden opportunity to convey our sincere gratitude to Asst. Prof. Debraj Das who helped us in carrying our project on COMPUTATIONAL FLUID DYNAMICS and provided useful guidance without which it would be really tough to complete this project. He was there in each and every stage to assist and motivate us, so that we could come up with a good work, and due to his faith and trust upon us, we were able to do this project work. This project also made us to know about the difference scope of CFD and the different governing equations and its link with physical physical problems of solidification solidification and melting of a material. material. At last we would like to thank our Head of Department Prof. Tarun Kanti Jana who has given us this opportunity to work with Sir. Debraj Das, and get an experience of his expertise in CFD.
Abstract In this study, basically we have dealt with Solidification and Melting problem, which is a moving boundary problem in which we track the solid- liquid interface which moves mo ves with time. Natural Convection and Conduction are the mechanism behind the physics of these problems. We have solved Navier-strokes equation along with continuity and energy equation, both in solid and liquid region using structured grid. In order to make zero velocity condition in solid domain special care has been taken. We have used enthalpy method to track the solid-liquid interface with respect to time. A fully coupled implicit method is used to solve the momentum and energy equation. A diffusion phase change, isothermal with convection along with continuous casting problem are present in the present study, and is validated with analytical and numerical results available. First, the two and three dimensional diffusion problem has been solved followed by gallium melting and mushy zone problem. Lastly, application problem on continuous casting has been solved and verified.
Contents 1
2
Introduction
2
1.1 1.2 1.3 1.4
2 4 7 7
Mathematical Modeling and Finite Volume Method
8
2.1 2.2
8 9 9 9 10 11 11 11
2.3
3
4
5
Methods needed for solving phase change problems Literature Survey Objectives Thesis Organization Assumptions Governing Equations 2.2.1 Continuity Equation 2.2.2 Momentum Equation 2.2.3 Energy Equation Initial and Boundary Conditions 2.3.1 Initial Conditions 2.3.2 Boundary Conditions
Results and Discussion
12
3.1
Diffusion Problem (Isothermal Case) 3.1.1 Two Dimensional Problem 3.1.2 Three Dimensional Problem 3.2 Isothermal Phase Change With And Without Convection 3.2.1 Gallium Melting with Convection 3.2.2 Gallium Melting without Convection(Diffusion) 3.3 Mushy Zone Problem 3.3.1 Two Dimensional Problem 3.4 Practical Application(Continuous Casting Process)
12 12 15 16 16 19 21 21 25
Conclusions And Scope For Future Work
31
4.1 4.2
31 31
Conclusions Scope For Future Work
References
32
List of Figures 3.1
Square Cavity Problem without Convection
12
3.2
Square Cavity Problem without Convection
13
3.3
Temperature distribution for square cavity problem (Case-1) (a)10 sec; (b)20 sec; (c)30 sec.
3.4
Position of interface in two dimensional problem (Case 2) (a) 41 sec; (b) 61sec; (c) 81sec.
3.5
14
14
Position of Interface (a) t= 0.60 sec, (b) t = 0.75 sec at z = 2 plane, (c) temperature contours at t = 0.75 sec at z
3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14
Melting of Gallium Problem Streamlines for Gallium melting (a) 6 min, (b) 9 min. Temperature Contours for Gallium Melting (a) 6 min, (b) 9 min Interface position at different times Melting of Gallium Problem Streamlines for Gallium melting (a) 6 min, (b) 9 min. Temperature Contours for Gallium Melting (a) 6 min, (b) 9 min Mushy region two dimensional problem Vector plot and mushy region for ϵ=0.1 (a) t=100 sec , (b) t=600 sec , (c) t=1000 sec 3.15 Comparision of (a) u velocity at t = 500 sec, (b) solidus and liquidus line at t =1000 3.16 Temperature contours for two dimensional mushy region problem (a) t =600 sec, (b) t =1000 sec. 3.17 Solidification in Czochralski Model
15 17 17 18 18 19 20 20 22 23 23 24 26
3.18 Shows the temperature contours for steady conduction solution
26
3.19 Shows Contours of Static Temperature (Mushy Zone) in steady State
27
3.20 Shows the Static temperature contour in transient state
27
3.21 Contours for Stream function at t = 0.2 sec.
28
3.22 Contours for liquid fraction
28
3.23 Contours of temperature t = 5 sec
29
3.24 Stream function Contours at t=5 sec
29
3.25 Contours of liquid fraction at t= 5sec
30
Nomenclature 3
A
Porosity function for the momentum equations (kg/m s)
b
Constant
C
Constant
e
Total Enthalpy (J/kg)
g
Acceleration due to gravity (m/s )
k
Thermal conductivity (W/mK)
p
Pressure (N/m2)
t
Time (s)
2
u, v, w
Velocity Component in x, y, z directions, respectively
c p
Specific heat at constant pressure (J/kgK)
eT
Sensible Enthalpy (J/kg)
e L
Latent Enthalpy (J/kg)
f l
Liquid fraction
Greek Letters
µ
Molecular Viscosity (kg/m-s)
Γϕ
Diffusion Coefficient of the variable ϕ
λ
Latent heat (J/kg)
1
Chapter 1 1. Introduction Now-a-days, the phenomenon of solidification and melting is of great importance in basic manufacturing processes like casting, welding etc. They have a great impact on many industrial applications. In earlier days, only analytical solutions were available, which did not give a clear idea about the process. Moreover, some effects (like convection) were also neglected in those days. So, implementation of numerical techniques for this kind of problems gathers attention for both present and future research. Solidification and melting problems are phase change problems, in which a solidliquid interface is moving with time and it has to be observed and tracked. One extra condition is required for solving general governing equations of this kind of problems. This condition is called Stefan condition and has to be applied at the solid-liquid interface. The Navier-Stokes equations coupled with the energy equation are solved in the problem domain. The problems can be solved numerically using computational fluid dynamics (CFD). In the present study, finite volume method (FVM) is used with structured meshes which can be easily applied in any arbitrary geometry. Ansys Fluent has been used as a tool to implement CFD, in the following thesis.
1.1 Methods needed for solving phase change problems There are many methods for solving the solidification and melting problems. As interface moves with time, they are classified according to the choice of domain. 1. Variable domain method : Here ,in this method, the governing equations are solved separately in both domains. Here the domain changes with time because the interface moves with time. For this reason it is called variable domain method. The Stefan‟s condition is applied to track the interface. So, this method requires adaptive grid generation and we have to track the interface. Two separate sets of equation for solid an d liquid are required. 2. Fixed domain method : Here, the domain does not changes with time. Governing equations are to be solved in the domain. The main disadvantage regarding this method is that it sometimes breaks down when interface moves a distance larger than a space increment in a time step. However, it can be easily solved using variable domain method. Again for solving multidimensional problems, variable domain method is not applicable.
2
There are number of techniques available for solving these problems and these can be found in Hu et al. [1]. These techniques are described below. (i). Strong solution method: In this method , for getting the interface for the new time step , the Stefan,s condition is applied directly and then governing equations are solved. This method requires large computation time in 3-D problems. Some methods under this category are of fixed grid method, level set method etc.
In this method, we do not use Stefan‟s condition directly instead (ii). Weak solution method: we use Stefan condition is implicitly incorporated in a new form of equations . These are also called as latent heat evolution methods. The main advantages of this method is that the explicit attention on the moving boundary is not required . There are different types of methods under this category, some of them are described below. a. Apparent heat capacity method: In this method, latent heat is accounted for by increasing the heat capacity of the material in the phase change temperature range. This method was first invented by Hashemi and Sliepcevich [2]. They used a finite volume formulation and Crank Nicolson for time integration. Apparent heat capacity is defined in three distinct region (like solid, liquid and mushy) depending upon the temperature. Although this method is computationally simple, it did not perform well compared to the other methods. For solving governing equation in case of a pure metal, we have to assume an artificial phase change temperature range which is a big disadvantage.
This method was proposed by Poirier and Salcudean [3] b. Effective heat capacity method: to eliminate the disadvantage regarding apparent heat capacity. In this case, a temperature profile is assumed in the interior nodes rather than calculating apparent heat capacity directly. The finite volume integration method is used to calculate the Effective heat capacity . The disadvantage of this metod is that it is difficult to implement and costly. c. Enthalpy method: In this case, the energy equation is transformed in the form of enthalpy. This transformed form of the energy equation i.e. enthalpy form is to be solved directly. Temperature of each cell is calculated from the enthalpy. For solving the phase change problem, it is the most efficient method and latent heat is incorporated through latent heat enthalpy term. Actually, total enthalpy is divided into two p arts namely sensible and latent enthalpy. Governing equation can be solved for sensible enthalpy with latent heat source term. Here, the relationship between enthalpy and temperature is important.
Generally, two types of phase change occur in metals and alloys. 3
Isothermal phase change:
In this case, phase change occurs at a distinct temperature,
enthalpy change is a steep change at melting temperature. This happens in case of pure metal i.e., Tin, Gallium etc. Mushy region phase change:
In this case, phase change occurs over a temperature range i.e.,
enthalpy becomes a continuous function of temperature. These problems are referred to mushy region phase change problems. The relationship between enthalpy and temperature can be any type linear, exponential. Here only linear relationship is considered. Binary alloys and all mixtures follow this relationship.
1.2 Literature Survey Earlier work related to solidification and melting problems is based on diffusion problems only and convection effects were not so dominant. A brief review regarding the modeling of solidification and melting problems can be found in Basu et al. [4] and Hu et al. [1]. Basu et al. [4] have described different types of methods such as fixed domain method, variable domain method for solving solidification and melting problems. They have formulated the governing equations for convection-diffusion phase change problems (isothermal as well as mushy region phase change case). Hu et al. [1] have formulated the governing equations through streamfunction-vorticity formulation as well as primitive variable formulation. Lazaridis [5] solved multi-dimensional diffusion problems by directly applying Stefan condition coupled with the energy equation. They solved four kind of diffusion problems. The discretization scheme for the cells surrounding the interface is different from that for the interior cells. They used both explicit and implicit time integration scheme. Voller and Cross [6] solved moving boundary problems using enthalpy methods. They used finite difference scheme for spatial and for time discretization both implicit and explicit methods are used. They solved two region problem and two-dimensional problem and compared the result with analytical and numerical result. Voller [7] developed implicit enthalpy formulation for binary alloy solidification without taking convection into account and used node jumping scheme for tracking solid-liquid interface. Crowley [8] extended multidimensional Stefan problems and he solved solidification of a square cylinder of fluid using enthalpy method when surface temperature is lowered at a constant rate.
4
Gau and Viskanta [9] first took the natural convection phenomenon in solidification and melting problems. They conducted an experiment for studying the buoyancy-induced flow in the melt and its effect on the solid-liquid interface position and heat transfer rate during the process of melt-ing and solidification of a pure metal (Gallium) from a vertical wall. They compared the solution with Neumann problem and concluded that convection effect can be neglected during phase change problems. Morgan [10] solved phase change problems taking convection into account. He used ex-plicit finite element method to solve freezing problem in a thermal cavity.
The basic enthalpy formulation of the governing equation was done by Voller et al. [11]. The enthalpy formulation is a weak solution method. They divided total enthalpy into sensible and latent enthalpy. They derived an equation for sensible enthalpy, in which latent enthalpy appeared as a source term. They solved the equation for sensible enthalpy and from that they calculated temperature. They used FVM for discretization. They solved a problem considering the effect of natural convection on isothermal solidification in a square cavity. They used different technique to make the velocity in solid region zero. They used variable viscosity method, Darcy source based method and switch-off technique as techniques for making zero velocity in solid domain. This approach was called “enthalpy- porosity technique”. Voller and Prakash [12] modelled a methodology for mushy region phase change problem by taking convection into account. They used enthalpy-porosity tech-nique as mentioned in Voller et al. [11] for formulation of governing equations. In mushy region, fluid velocity is not zero and therefore mushy region contributes to some convection, they assumed that in mushy region flow occurs through a porous media. They defined „permeability‟ to model the flow and they took same governing equation which relates fluid velocit y and pressure, derived from the Darcy law.
− − (
=
)
(1.1)
where λ is the permeability of the porous medium. Voller et al. [11] neglected convective latent enthalpy source term for isothermal phase change case. Voller and Prakash [12] did not neglect the convective term of latent enthalpy source term as it is not zero in case of mushy region phase change problem. They derived general formulae for both temporal and convective latent enthalpy source term. Brent et al. [13] applied the formulation proposed by Voller and Prakash [12], to the problem of the melting of Gallium in a rectangular cavity. They considered isothermal case and convection was taken into account. They plotted isotherms and streamlines at different times and compared their results with the experimental results obtained by Gau and Viskanta [9]. Wolff et 5
al. [14] solved problem regarding the solidification of Tin in a square cavity by using numerical as well as experimentally. The two sides of the cavity were at a fixed temperature and remaining two were insulated. At last they compared the numerical result with experimental result. For numerical technique they used enthalpy method. Rady and Mohanty [15] used enthalpy-porosity technique to solve melting of Gallium . They validated their result with Wolf et al. [14]. They plotted isotherms and streamlines at different times in case of melting of Gallium problem. They also plotted the interface position at different times. Stella and Giangi [16] studied the melting of pure Gallium in a bi-dimensional rectangular cavity with aspect ratio 1.4. They plotted solid-liquid interface and streamlines at different times and shown a multi-cellular flow structure built in the process of melting. Redy et al. [17] studied about the effects of liquid superheat during solidification of pure metals. They also used the enthalpy-porosity technique. They obtained steady state very early for higher Rayleigh numbers. They plotted Nusselt number variations and temperature profiles for different Rayleigh numbers. Ghasemi and Molki [18] studied isothermal melting of a pure metal enclosed in a square cavity having Drichlet boundary conditions in each side. They continued their 8
7
computations for Rayleigh number 0 to 10 and Archimedes number 0 − 10 . They plotted liquid fraction variation with time, falling velocity of solid phase and shape of the solid-liquid interface. They found that for low Rayleigh and Archimedes number, both melting rate and solid velocity are low and melting is almost symmetrical. Melting rate enhances with the higher value of Rayleigh and Archimedes number.
Gong and Mujumder [19] studied melting of a pure phase change material in a rectangular con-tainer heated from below. They used Streamline Upwind/Petrov Galerkin finite element method in combination with fixed grid primitive variable method. Flow patterns for different Rayleigh numbers were used. They also studied the instability of free convection flow at higher Rayleigh numbers. Bertrand et al. [20] reviewed the methods to solve the solidification problems and compared the results. They gave the results for high as well as low Prandtl number fluids. Hwang et al. [21] considered the effect of density variation with phase change when tin solidifies in a square cavity. They used multi-domain method to cope up with abnormal variations of front position due to shrinkage.
6
1.3 Objectives The objective is to solve Three and Two Dimensional Solidification and Melting Problem using Ansys Fluent 16.2 for both Isothermal and mushy region phase change and validate the simulation results with numerical, experimental and analytical solutions available in the literature.
1.4 Thesis Organization A brief introduction along with literature review is presented in chapter 1. Mathematical modeling, Governing equations and initial and boundary conditions are described in chapter 2. Problem solving using Ansys Fluent on Solidification and Melting is shown in chapter 3. At last Conclusion and scope for future works are listed.
7
Chapter 2 2. Mathematical Modeling and Finite Volume Method Now – a-days, Navier-Stokes equations and energy equation in solidification and melting problems are solved using the fixed domain method. Because of versatility of Fixed domain enthalpy method ,it can be used for both isothermal and mushy region phase change problems. In this case, as the position of the interface is obtained as part of the solution, explicit information about the interface is not required. While solving Navier-Stokes equation in the solid domain, attention must be taken to make zero velocity condition in that domain. Therefore, the fixed domain enthalpy method demands some techniques to do that, which is described in the next. As convective effect is not neglected so the Navier-Stokes equations and the energy equation are coupled in these problems . In the present formulation, the governing equations have been considered in Cartesian coordinates system.
2.1 Assumptions 1. The flow is considered to be incompressible, Newtonian and laminar. 2. Properties like thermal conductivity, specific heat are assumed to vary linearly with liquid fraction. 3. The density variation due to phase change is neglected for closed domain problems (like square cavity problem). The density variation due to temperature in the liquid domain is incorporated through Boussinesq approximation. Variable density formulation is to be used in case of external flow. However, the variable density formulation cannot handle shrinkage effect during solidification. This needs some special treatment [21]. 4. Species transport equation is not solved, so solute buo yancy is not included. Only thermal buoyancy is considered in the present study. 5. Viscous dissipation effect is neglected.
8
2.2 Governing Equations The governing equations are as follows based on above assumption are written below.
2.2.1 Continuity Equation
.(
)=0
( . )
2.2.2 Momentum Equation The Navier-Strokes equation (in vector form) for laminar, incompressible flow of Newtonian fluid can be written as follows
∇ −∇ ∇ μ∇ ρ (
)
+ .
=
p+
.
+
( . )
To make velocities equal to zero in the solid domain, a large negative source term is added to the above equation. The source term becomes zero when it is liquid domain. So, the equation then becomes
∂ ∂ ∇ −∇ ∇ μ∇ ρ (
)
t + .
=
p+
.
+
+A
( . )
The second source term takes very high value for making the velocities very close to zero in the solid domain and in the liquid domain, it is simply zero. The equation for A is [13]
− −
( . )
( . )
)2
(1
=
3
+
Where f l is the liquid fraction ,which is defined as the ratio of volume of the liquid present in any particular cell to the total volume of the cell.
=
„C‟ and „b‟ arre prescribed constants. The equation for „A‟ makes the momentum equation to follow the Carman-Kozeny equation in the mushy region. In mushy region both solid and liquid phase are present and therefore f l always takes the value like 0< f f <1. Therefore, it is assumed that the flow occurs through a porous media [13]. Both „A‟ and „C‟ have dimension to satisfy
) and „b‟ is assigned a value 0.001 to avoid division by zero. The value of both „C‟ and „b‟ makes very high value of „A‟ in the solid 6
momentum equations and „C‟ has a very high value (10
3
region and fluid velocity becomes close to zero as f l =0. But in liquid region f l =1 and therefore the entire source term becomes zero.
9
2.2.3 Energy Equation The general form of energy equation after neglecting viscous dissipation term is Eqn 2.6 The direct form of enthalpy equation is used in the present work. Enthalpy is split into two parts i.e.,
=
+
( . )
Where eT is sensible enthalpy and eL is latent enthalpy per unit mass.
=
( . )
eL=0 in the solid region, eL =L in the liquid region and eL varies between 0 and L for the cells undergoing phase change. Substituting e in the energy equation, ( + ) + . + = . ( . )
∇ ∇ ∇ ∇ ∇ ∇ − −∇
After simplifying the equation becomes
(
)
+ .
(
= .
)
.
( .
)
Now, as the right hand side of the above equation is in terms of temperature, its contribution to the diagonal term coefficient is zero. So the equation becomes similar to one of pure convective equation with the source term, which is less stable during numerical solution. To improve convergence rate and stability of the above equation, temperature in diffusion is replaced b y eT. Then the equation becomes
∇ ∇ ∇ – −∇ (
)
+ .
=
.
.
( .
Latent enthalpy term in the right hand side can be considered as source terms. Detailed updation procedure can be found in [12].
10
)
2.3 Initial and Boundary Conditions 2.3.1 Initial Conditions Initial conditions combined with the boundary conditions will determine whether the given problem is a phase change problem or not. As phase change problems are unsteady problems, the initial conditions play a major role in the solution. For solidification, initially some part of the domain has to be liquid. According to the types of problems parameters are to be initialized. While solving energy equation, all initial and boundary conditions have to be in terms of sensible enthalpy.
2.3.2 Boundary Conditions Boundary conditions needed for the solidification and melting problems can be Dirichlet, Neumann or Robin. The boundary condition is to be implemented as follows. For Neumann boundary condition is implemented as,
− =
( .
)
For solid walls no-slip boundary conditions are used. For pressure, the homogeneous Neumann boundary condition is used for velocity specified boundaries. For other boundaries, appropriate boundary conditions should be specified to the physics of the problem.
11
Chapter 3 3. Results and Discussion In this chapter, problems related to solidification and melting are discussed and verified with the numerical and analytical solution. Firstly, diffusion phase change problems without convection effect (Isothermal case only) and then both isothermal and mushy region convection-diffusion phase change problems are discussed. For solving, phase change problems with convection, a suitable source term are added in the momentum equations to get zero velocity in the solid domain. The sensible enthalpy form of the energy equation is solved with the appropriate source terms. Lastly, a problem on continuous casting has been also discussed in the given thesis, which portrays the application part of this solidification and melting problem.
3.1 Diffusion Problem We have solved two benchmark problems on solidification and melting to corroborate the simulation done through Ansys Fluent, with the numerical and Analytical solution from the literature. The first question being a 2D problem and the other next one is the 3D problem. Efforts have been made to simulate and bring the result in accordance with the literature. The properties taken in such a way that they matches with some non-dimensional quantities (e.g. Stefan number or some parameter defined in the literature).
3.1.1 Two- dimensional problem This is a two dimensional problem in which solid is melted in a square cavity having same wall temperature. There are two cases for solving this problem. o
Case 1:-Cavity wall temperature is 1 C.
Figure 3.1:- Square Cavity Problem without Convection 12
Case 2:- Cavity wall temperature is 0.5o C
Figure 3.2:- Square Cavity Problem without Convection
For solving this problem we have to take a square cavity having 1x1 dimension. Here four interfaces are formed i.e. four solid walls and these are joined to form a single interface. In both 0 cases, initially the solid is kept inside the cavity at its melting point i.e. Tm= 0 C. After that the temperature of all the boundary walls are increased suddenly. Different physical properties taken for this square cavity problem are as follows:-
Table 3.1:-Physical properties taken for square cavity problem
k
C p
ρ
λ
µ
β
(W/m-K)
(J/kg-K)
(kgm )
(J/kg)
(kg/m-s)
(1/K)
Solid
1.0
100.0
1.0
0.0
-
-
Liquid
1.0
100.0
1.0
1000.0
0.1
0.01
-3
Now for computation purpose, we have to choose a 81 x 81 grid mesh. And the above given 4 physical properties chosen such that these match with non-dimensional parameters like Ra = 10 , St=0.1 ,Pr =10.
13
(a)
(b)
(c)
(d)
Figure 3.3:- Temperature distribution for square cavity problem (Case-1)(a)10 sec; (b)20 sec; (c)30 sec.
(d)40 sec
(a)
(b)
(c)
Figure 3.4:- Position of interface in two dimensional problem (Case 2) (a) 41 sec; (b) 61sec; (c) 81sec
14
The solution graph shows interface position at different time inside the cavity. It also shows that Temperature/ Contour lines are axis-symmetry as the boundary conditions of each faces is same and the geometry is symmetric. Here the intensity of convection is less as Raleigh number is less 4, I.e Ra =10 and as we know that as Rayleigh number increases, the tendency of pushing the solid increases. And here, Rayleigh number is less so liquid does not have enough potential to push the solid. So in this case conduction i.e diffusion of heat dominants over natural convection as Rayleigh number is less. so , convection is neglected in this problem .
3.1.2 Three-dimensional problem As we have discussed about the 2-dimensional solidification melting problem previously, now we are getting interest to study about the simulation of the solidification-melting problem of a 3-dimensional cavity having dimensions (4x4x4). This benchmark problem is taken from the thesis 0 Lazaridis [5] .Initially in this problem the liquid metal having melting point 0 C is kept in the 0 cavity. Suddenly, the temperature of the of the left and bottom wall is reduced to -3.24 C that means Dirichilet condition is applied and the other four boundaries having Neumann boundary conditions which means they are adiabatic in nature. Now we will track the position of the liquid solid interface at a constant plane that is Z=2 at different times and will compare with the thesis of Lazaridis [5]. During the simulation a time step of .01 is chosen. The properties of the metal are taken by using some non-dimensional number according to Lazaridis [5]. The properties are given in the Table 3.2. Table 3.2:- Properties for 3D problem
Solid Liquid
(a)
k (W/m-K)
C p (J/kg-K)
ρ -3 (Kgm )
λ (J/kg)
1 1
1 1
1 1
0 5
(b)
(c)
Figure 3.5-: Position of Interface (a) t= 0.60 sec, (b) t = 0.75 sec at z = 2 plane, (c) temperature contours at t = 0.75 sec at z = 2 plane.
15
Figure 3.5(a) and (b) is the comparison of the interface of our simulation with the thesis of Lazaridis [5] at different timings. Figure 3.5(c) is showing the temperature contours at Z=2 and t=0.75 sec. Here the effects of conduction are much higher than the effect of convection so we consider only the effect of conduction. Here we observe that the wall having lower temperature converted to solid very quickly then the position interface progress in such a way showed the Figure 3.5. The interface is looking like a parabola.
Table 3.3-:Physical Properties of Gallium
Solid Liquid
k (W/m-K)
C p (J/kg-K)
ρ -3 (Kgm )
λ (J/kg)
µ (kg/m-s)
β (1/K)
32.0 32.0
381.5 381.5
6095.0 6095.0
0.0 80160.0
-3 1.81x10
-4 1.2x10
3.2 Isothermal phase change with and without convection The phase change problems solved till now, only deals with diffusion. Since, the research says that convection effect cannot be neglected [9], therefore the problems now dealt with are solved with convection effect considered and variation obtained is studied and analyzed.
3.2.1 Gallium Melting with Convection As we know Gallium melting, being a benchmark problem for isothermal phase change problems, convection effect was considered and the simulation was validated using Ansys Fluent. A problem from Brent et al. [13] is taken to evaluate the simulation obtained. The problem is defines as stated below: A rectangular cavity of 0.0889 m in length and 0.0635 m in height is taken in which pure solid °
°
Gallium is initially kept at 28.3 C. Suddenly, left wall temperature is increased to 38 C which is higher than the melting point temperature (T m= 29°C) of Gallium and the right wall is kept at the initial temperature of the solid Gallium. Other two boundaries are insulated. A 42 x 32 grid is chosen for the simulation. Figure 3.6 shows the computation domain with the necessary boundary conditions. The physical properties are taken from Brent et al [13]. The properties are shown in the 6 table 3.3. For the present simulation, A=10 and b=0.001 are taken. A time step of 0.01 is used.
16
Figure 3.6: Melting of Gallium Problem
(a)
(b)
Figure 3.7: Streamlines for Gallium melting (a) 6 min, (b) 9 min Figure 3.7 shows streamlines at different times. Melting of Solid Gallium takes place due to the heated wall and a solid-liquid interface moving in the right hand side of the cavity and therefore, density of the liquid Gallium changes with temperature in the area adjacent to the heated wall. The liquid Gallium rises up having less density and heavier liquid stays at the bottom. As a result, a convection current due to density difference in liquid Gallium is set up inside the cavity, which is called as natural circulation, and this enhances the melting. Figure 3.8 shows temperature contours at different times. Initially, the contours are straight which indicates that heat transfer occurs mainly due to conduction. But as the time progresses, convection phenomena becomes dominant and a slight curvature in the contour plot is observed.
17
Figure 3.9 shows interface comparison with Brent et al. [13] at different times. Front position comparison obtained is satisfactory. To plot the interface position, f l =0.5 contour is used. At any particular time, the interface divides the cavity area into two distinct phases. With the passage of time, front moves rapidly in the top but the movement is very slow in the bottom. This signifies that the convection is more prevalent in the upper portion. Hot liquid impinges on the solid Gallium in the top of the cavity and therefore melting rate is more in the top as compared to
(a)
(b)
Figure 3.8:-Temperature Contours for Gallium Melting (a) 6 min, (b) 9 min
Figure 3.9: Interface position at different times
18
the bottom portion. However, melting of Gallium is a controversial problem in the literature [26]. There is a controversy between multi-cellular and mono-cellular liquid flow of Gallium, first pointed out by Dantzig [27]. In the present study, mono-cellular liquid flow is assumed. A great effort has been made by Hannoun et al. [26] to solve the controversy. However, the interface position does not match well with the literature. It overestimates the result given in Brent et al, the reason could be time step used and discretization techniques. Due to controversial effect discussed in the previous paragraph, may also play an important role in the estimation of the interface position at different times.
3.2.2 Gallium Melting Without Convection Gallium melting as already stated, being a benchmark problem for isothermal phase change problems, in which now the convection effect earlier considered and was neglected and the suitable data value are obtained and studied. Again the same problem from Brent et al. [13] is dealt with convection effect neglected. For the convenience to reader again the problem is defined in the same manner as it was done in the previous problem.
Figure 3.10:- Melting of Gallium Problem Similar post processing is done as in Gallium melting with conve ction problem, with the aim to correlate both the problems deduce the dominating factor prevalent in the problems.
19
(a)
(b)
Figure 3.11:Interface for Gallium melting (a) 6 min, (b) 9 min Figure 3.11 shows interface at different times. Melting of Solid Gallium takes place due to the heated wall and a solid-liquid interface moves in the right hand side of the cavity. Density of the liquid Gallium changes with temperature in the area adjacent to the heated wall. Since the convection effect is neglected no circulation phenomenon is developed. Figure 3.12 shows temperature contours at different times. All through the simulation period, contours are straight which indicates that heat transfer occurs mainly due to conduction. Since, only convection predominates only straight contours are obs erved, with negligible streamlines.
(a)
(b)
Figure 3.12: Temperature Contours for Gallium Melting (a) 6 min, (b) 9 min
20
3.3 Mushy region problem Till now, we have dealt with isothermal phase change problems with and without natural convection. In case of isothermal phase change, phase change occurs at distinct temperature. But, in case of mushy region, phase change occurs over a temperature range. Therefore, latent heat is dependent upon the temperature. Several relationships between latent heat and temperature can be possible. A linear relationship is mostly used in literature due to its simplicity. In the present case, both linear and non linear relationships are taken. In isothermal phase change, we have neglected the convective term of the latent heat, but it is included in mushy region [12]. Both two and three dimensional problem are presented here. The linear relationship can be obtained as follows:
=
+
(3.1)
Where a and c are constants. Now, the value of the constants can be determined by applying suitable conditions. At T=Ts, e L =0 and at T=TL, e L=λ. By applying these condition the following relationship is obtained
− − =(
(
) ;
(3.2)
A more general relation is
λ − − =(
) ;
(3.3)
Where n is the index value. Generally, 2 ≤ n ≤ 5 is accepted. When n=1, the relationship becomes linear. Otherwise, it is non-linear.
3.3.1 Two-dimensional problem A two dimensional mushy region phase change problem is taken from Vollar and Prakash [12]. However, it has been reported in [24] that for high Prandtl number ( Pr=1000) liquid, the requirement of computational time is more. Therefore, to reduce the computation time, low Prandtl number ( Pr = 10) is taken by Debraj Das [24]. We have also taken Pr =10 in the simulation. The computational domain is shown in Fig. 3.13. A 40 X 40 uniform hexahedral grid is taken in 1 x 1 domain. Initially, cavity is filled with liquid having initial temperature 0.5◦C. Suddenly, left wall temperature is reduced to -0.5◦C and right wall is kept as the initial temperature 0.5◦C. Suddenly, left wall temperature is reduced to 21
0.5◦C and right wall is kept as the initial temperature. Other two boundaries are kept insulated. The problem is solved with different values of ϵ .The value of the liquid properties shown in Table 3.4 4 are calculated keeping the non- dimensional number same ( Ra=10 , Pr =10, St =5.0) [24]. In the 3 present computation, the value of constant A and b is taken as 10 and 0.01. Figure 3.14 shows
Figure 3.13:-Mushy region two dimensional problem Table 3.4:-Physical Properties taken for two dimensional mushy region problem
Solid Liquid
k (W/m-K)
C p (J/kg-K)
Ρ (Kgm )
Λ (J/kg)
µ (kg/m-s)
β (1/K)
T f (K)
0.001 0.001
1.0 1.0
1.0 1.0
0.0 5.0
0.01
0.01
0.0
Vector plot along with the position of mushy region at different times. For simulation validation the value of ϵ is taken as 0.1. The mushy zone is the region bounded by the solidus and liquidus lines. Therefore, region bounded by T L = ϵ and T S = - ϵ is the mushy zone as ϵ = (T L- T S )/2. It is seen from the vector plot that fluid velocity is not zero inside the mushy zone, so the mushy zone contributes some convection which is expected. As time progresses mushy zone moves towards the right wall and the solidification rate is increased. However, in solidification there can exist a steady state [24]. In this problem, it is seen that mushy zone does not move very much with time after t = 1000 sec. So, the problem reaches steady state. In the present section, unsteady results
22
(a)
(b)
(c)
Figure 3.14:- Vector plot and mushy region for ϵ =0.1 (a) t=100 sec , (b) t=600 sec , (c) t=1000 sec are presented only. Comparison of u velocity at different horizontal sections at t = 5 00 sec is shown in Figure 3.14(a). Figure 3.14(b) shows the comparison of solidus and liquidus line at t =1000sec. The comparison is good. The velocity variations are sinusoidal in nature. Small Variation in results is found due to difference in grid. Figure 3.16 shows the temperature contours at different times. Near to the left wall contours are straight which indicates that the heat transfer mainly occurs due to conduction only. Figure 3.15(b) shows the effect of half mushy range (ϵ ) on the width of the mushy region at t = 1000 sec.
(a) (b) Figure 3.15:- Comparision of (a) u velocity at t = 500 sec, (b) solidus and liquidus line at t =1000
23
(a)
(b)
Figure 3.16:-Temperature contours for two dimensional mushy region problem (a) t =600 sec, (b) t =1000 sec The problem has been solved by assuming a linear relationship between latent enthalpy and temperature.
24
3.4 Practical Application (Continuous casting process) Continuous casting process is now-days, a backbone manufacturing process for almost all the manufacturing firms, especially the steel plants. With advancement in technology, we have seen the use of enhanced and improved manufacturing methods, to eliminate earlier defects like porosity, blow holes, air inclusion, slag entrapment and many more. Here, therefore we have implemented Ansys Fluent to simulate the real life problem virtually, to effectively study the outcomes in order to draw conclusions and eradicate the problems faced in manufacturing industries. Here, we are solving a benchmark problem on Solidification and melting commonly known as “Solidification in Czochralski Model”, which takes place in Continuous Casting Process and is directly taken from Ansys Fluent tutorial, defined as a solidification problem which involve fluid flow and heat transfer problem. A 2D axis symmetric bowl has been considered as geometry as shown in the Figure 3.17 which contain the liquid, the boundary conditions have been mentioned on the figure itself. The bottom and the sides of the bowl are heated above the liquidus temperature, as it is the free surface of the liquid. The liquid is solidified by heat loss from the crystal and the solid is pulled out of the domain at a rate of 0.001 m/s and a temperature of 500 K. A steady injection of liquid is maintained at the bottom of the bowl with a velocity of 1.01 x 10
-3
m/s and a temperature of 1300 K. Material properties. Initially, steady solution is computed and then the fluid flow is enabled to investigate the effect of natural and Marangoni Convection in a transient fashion. As in continuous casting the material is pulled out continuously so, pull velocity is enabled. A plot of temperature contours is shown in Figure 3.18, in accordance with the literature present in Ansys Fluent Tutorial [23]. Contours of temperature (Mushy zone) is shown in Figure 3.19-: in steady state solution,
25
Figure 3.17: Solidification in Czochralski Model Case 1: Steady State Solution
In steady state solution we would specify the type of discretization form, and would restrict the calculation of flow and swirl velocity equation on temporary basis, in order to calculate conduction only.
Figure 3.18:-Shows the temperature contours for steady conduction solution.
26
Figure 3.19:- Shows Contours of Static Temperature (Mushy Zone) in steady State. Case 2: Transient State Solution
The previous steady state solution was taken as the initial condition for transient state solution. In this flow and swirl velocity equation is being calculated. Now, for transient state solution numerous plots have been plotted in various parameters to validate the problem with the thesis [23].
Figure 3.20:- Shows the Static temperature contour in transient state
27
Figure 3.21:- Contours for Stream function at t = 0 .2 sec
In the Figure 3.21 , the contour obtained for Stream function is due to Natural convection and Marangoni convection on the free surface.
Figure 3.22-:Contours for liquid fraction
Figure 3.22, here the current position of melt front is being displayed here. The Mush y zone divides the liquid and solid regions roughly into half. Further simulation is done for 48 more number of time steps in all total of 50 time steps. And the plots obtained us somewhat like this. In Continuous Casting, it is very important to know when to pull the material out. If it is pulled earlier it won‟t get solidified (still in mushy zone), and if pulled late it will be solidified and won‟t be pulled out in a desired shape. By proper study of the solidus and liquidius temperature contour, the optimal pull rate can be easily calculated. 28
Figure 3.23:- Contours of temperature t = 5 sec
The Figure 3.23, shows that the temperature contour is fairly uniform along the melt zone and solid material, and due to recirculating liquid the distortion obtained in the temperature is clearly visible. After 5sec the flow is mainly dominated by Natural convection and Marangoni stress as seen from Figure 3.24.
Figure 3.24: Stream function Contours at t=5 sec
29
Figure 3.25:- Contours of liquid fraction at t= 5sec
30
4. Conclusions and scope for future work 4.1 Conclusion The following study entrails the simulation done through Ansys Fluent, on few benchmark solidification and melting problems which include 2D and 3D isothermal problem, Gallium melting with and without convection, a Mushy region and problem and in the last a application problem on Continuous casting problem named as “Solidification in Czochralski Model”, and has been validated with the numerical and analytical data present in the literature. In the first two problems for isothermal case (diffusion problem) has been solved, followed by mushy region problem. In Gallium melting with and without convection has been tried to solve to show that convection plays a vital role in enhancing the rate of melting. Many Plots on temperature, streamlines and interface has been plot to show the variation obtained. Natural Convection pattern has been shown using the streamline plot. In the last problem, on continuous casting, both steady state and transient state solution has been shown.
4.2 Scope for Future Work We have neglected Density Change in the present thesis. With rapid void and Shrinkage formation because of Density change the flow field changes to shrinkage-induced flow. In the present study species equation has not been solved. So, one could add with Navier-Strokes and Energy equation for determining the composition in each phase in the multi-component alloy solidification process. Only thermal buoyancy has been considered. However, due to nonuniform cooling solutal buoyancy effect may arise. The flow has been considered as laminar. Models related to turbulent flow can also be added to phase change problems. The Ansys Simulation results has been verified with the literature data.
31
5. References [1] H. Hu and S. A. Argyropoulos, “Mathematical modelling of solidification and melting,” Modelling Simul. Mater. Sc. Engineering, vol. 4, pp. 371 – 396, 1996. [2] H. T. Hashemi and C. M. Sliepcevich, “A numerical method for solving two-dimensional problems of heat conduction with change of phase,” Chemical Engineering Progress Symposium, vol. 63, no. 79, pp. 34 – 41, 1967. [3] D. Poirier and M. Salcudean, “On numerical methods used in mathematical modeling of phasechange in liquid metals,” ASME Journal Heat Transfer, vol. 110, pp. 562 – 570, 1988. [4] B. Basu and A. W. Date, “Numerical modelling of melting and solidification problems-a review,” Sadhana, Part 3, vol. 13, pp. 169– 213, 1988. [5] A. Lazaridis, “A nemerical solution of the multidimensional solidification (or melting) problem,” Int. J. Heat Mass Transfer, vol. 13, pp. 1459– 1477, 1970. [6] V. R. Voller and M. Cross, “Accurate solutions of moving boundary problems using enthalpy method,” Int. J. Heat Mass Transfer, vol. 24, pp. 545– 556, 1981. [7] V. R. Voller, “An implicit enthalpy solution for phase ch ange problems:with application to a binary alloy solidification,” App. Math. Modelling, vol. 11, pp. 110– 116, 1987. [8] A. R. Crowley, “Numerical solution of stefan problems,” Int. J. Heat Mass Transfer, vol. 21, pp. 215 – 219, 1978. [9] C. Gau and R. Viskanta, “Melting and solidification of a pure metal on a vertical wall,” Journalof Heat Transfer, vol. 13, pp. 297 – 318, 1988. [10] K. Morgan, “A numerical analysis of freezing and melting with convection,” Computer Methods in Applied Mechanics and Engineering, vol. 28, pp. 275 – 284, 1981. [11] V. R. Voller, M. Cross, and N. C. Markatos, “An enthalpy method for convection/diffusion phase change,” Int. J. Num. Methods Engg, vol. 24, pp. 271 – 284, 1987. [12] V. R. Voller and C. Prakash, “A fixed grid numerical modelling methodology for convectiondiffusion mushy region phase change problems,” Int. J. Heat Mass Transfer, vol. 30, no. No.8,pp. 1709 – 1719, 1987. [13] A. D. Brent, V. R. Voller, and K. J. Reid, “Enthalpy-porosity technique for modeling convection-diffusion phase change: Application to the melting of a pure metal,” Numerical Heat Transfer, vol. 13, pp. 297 – 318, 1988. [14] F. Wolff and R. Viskanta, “Solidification of pu re metal at a vertical wall in the presence of liquid superheat,” Int. J. Heat Mass Transfer, vol. 31, pp. 1735– 1744, 1988. [15] M. A. Rady and A. K. Mohanty, “Natural convection during melting and solidification of pure metals in a cavity,” Numerical Heat Transfer P art A, vol. 29, pp. 49 – 63, 1996. [16] F. Stella and M. Giangi, “Melting of a pure metal on a vertical wall: Numerical simulation,” Numerical Heat Transfer Part A, vol. 38, pp. 193 – 208, 2000. [17] A. M. Rady, V. V. Satyamurty, and A. K. Mohanty, “Effects of liquid superheat during solidification of pure metals in a square cavity,” Heat and Mass Transfer, vol. 32, pp. 499– 509, 32
1997. [18] B. Ghasemi and M. Molki, “Melting of unfixed solids in square cavities,” International Journal of Heat and Fluid Flow, vol. 20, pp. 446 – 452, 1999. [19] Z. Gong and A. S. Mujumdar, “Flow and heat transfer in convection dominated melting in a rectangular cavity heated from below,” Numerical Heat Transfer Part B, vol. 23, pp. 21– 41, 1993. [20] O. Bertrand, B. Binet, H. Combeau, S. Coutturier, Y. Delannoy, D. Gobin, M. Lacroix, P. L. Quere, M. Medale, J. Mencinger, H. Sadat, and G. Vieira, “Melting driven by natural convection a comparison exercise:first results,” Int. J. Them. Sci, vol. 38, pp. 5– 26, 1999. [21] K. Hwang, “Effects of density change and natural convection on the solidification process of a pure metal,” J. Materials Processing Technology, vol. 71, pp. 466– 476, 1997. [22] A. Dalal, V. Eswaran, and G. Biswas, “A finite-volume method for navier-stokes equation on unstructured meshes,” Numer ical Heat Transfer, Part B, vol. 52, pp. 238 – 249, 2008. [23] Ansys Tutorial 2008 Full Guide, Melting and Solidification Problem, Chapter 24, “Solidification in Czochralski Model” [24] Debraj Das, “Numerical simulation on Solidification and Melting Problems on Unstructured Grid.”.
33
