- Email: [email protected]
Numerical study of melting a rod by a periodically moving local heat source
International Journal of Thermal Sciences 97 (2015) 1e8
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Numerical study of melting a rod by a periodically moving local heat source Vadim Mizonov a, *, Nicolay Yelin b, 1 a b
Department of Applied Mathematics, Ivanovo State Power Engineering University, Rabfakovskaya 34, 153003, Ivanovo, Russia Department of Hydraulics, Heat and Mass Transfers, Ivanovo State Polytechnic University, 8 Marta 20, 153037, Ivanovo, Russia
a r t i c l e i n f o
a b s t r a c t
Article history: Received 21 July 2014 Received in revised form 14 June 2015 Accepted 15 June 2015 Available online xxx
The objective of the study is developing a simple yet informative mathematical model that describes the kinetics of melting a rod under the action of a localized periodically moving heat source. For this purpose, a cell model is used with the heat conduction matrix that takes into account different properties of liquid and solid phases of the rod material. Zones of the rod that are outside of the local heat source action have heat exchange with the outside environment. It is shown that the melting kinetics strongly depends on the program of heat source motion along the rod and on its residence time at each of its positions. The optimal program of heat source motion that allows melting the whole rod over the shortest possible period of time is found. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Heat conduction Melting Solidification Local heat source Heat transfer Optimal program
1. Introduction The problem in question arises in some technologies when there is a need to heat up and melt an object by a heat source whose dimensions are much smaller than those of the object to be melted. The simplest example is defrosting a finite fragment of a water pipe full of frozen water by a burner. Sometimes (for instance, at a low outside temperature), this cannot be done if the heat source is only applied to one particular point of the fragment, and it is only its motion along the pipe that can solve the problem. In this case, the complete defrosting time strongly depends on the program of the burned motion over the pipe length e which can be the objective function of investigation. Analogous problems can arise in other technologies, for example, in hand treatment of roof materials by a gas burner, in soldering, and in additive technologies when laser beam treatment of thin layer of fine powder occurs. This problem is related to the non-linear heat conduction problems with a localized moving heat source. Its non-linearity is
* Corresponding author. Tel.: þ7 910 9948858; fax: þ7 4932385701. E-mail addresses: [email protected] (V. Mizonov), [email protected] (N. Yelin). 1 Tel.: þ7 910 9810020. http://dx.doi.org/10.1016/j.ijthermalsci.2015.06.005 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.
conditioned by the phase transformation of melting/solidification type and various thermal effects connected with it. Different aspects of the problem solution can be found in literature. However, most of them are connected with the analytical or semi-analytical approach to solving this problem e which inevitably implies using far-going assumptions. A solution to the problem of heat conduction in a rectangular plate exposed to a moving heat source was presented in Ref. [1]. The heat source moved along an elliptical trajectory that always remained within the boundaries of the plate area. The exact solution to the problem in an analytical form was obtained by applying the Green's function method. Exemplary results of numerical calculations to determine the temperature distribution in the plate were presented. However, first of all, this was the only trajectory of the source motion and, second, the phase transformation was not taken into account. The problem of the object being heated by a moving source with application to welding was also examined in Ref. [2] where the meshless local Petrov-Galerkin method was developed. The problem of bead-on-plate welding (a moving heat source problem), and the fundamental properties of the method were investigated to verify the applicability of the proposed method. The investigation was mainly oriented to the influence of the number of nodal points on the accuracy of the solution. The issue of the influence of the program of the source motion on the
2
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
process properties was beyond this particular study. The study presented in Ref. [3] describes the analytical and numerical solution of the heat conduction equation for a localized moving heat source of any type to be used in laser material processing, such as welding, layered manufacturing and laser alloying. The analytical solution for a uniform heat source was derived from the solution of an instantaneous point heat source. The result was evaluated numerically and was compared to existing solutions for the moving point source and a semi-ellipsoidal source. Next, the result was used to demonstrate how such model can be used to study the effect of the heat source geometry. To investigate the effect of the temperature dependence of the material parameters, in particular the latent heat of fusion, a finite difference model was implemented. In spite of the fact that many components of the process were taken into account the influence of the program of the source motion remained beyond the scope of attention of the study. An analytical approach of transient heat conduction in two-layered material, of finite depth, with an imperfect thermal contact, subjected to a moving gaussian laser beam was developed in Ref. [4]. The method consists in deriving the solution of the homogeneous part of the heat equation by using the well known separation of variables method and expressing the source term in series form. This model has been successfully applied on a practical system; laser cladding of electronic copper tracks on alumina substrates. This analytical model can also be used for estimation of the thermal contact resistance between the layers. However, the model does not take into account the phase transformation. An analytical method of computation of temperature field in a half-infinite body caused by heat source with changeable direction of motion was presented in Ref. [5]. Analytical temperature field was approximated by straight segments for volumetric heat source with a trajectory considering temperature changes caused by next transitions (increase in temperature connected to action of the heat source and self-cooling of areas heated-up earlier). In that instance computations were carried out for cuboidal elements made of steel for various heat trajectories. However, the phase transformation in the body was not taken into account too. Later on, in the paper [6], the melting process was introduced in the model but only as a reason of the increase of heat transfer from the heat source to the body. Melting of the body itself was not examined. A model that describes the transient heating of a thin wire causing the tip to melt, roll-up of the molten mass into a ball due to surface tension forces, and the subsequent solidification of the molten material due to conduction up the wire and convection and radiation from the surface, was proposed in Ref. [7]. The wire was assumed to be heated at its lower tip to a temperature beyond the melting temperature of the wire material by heat flux from an electrical discharge. However, the objective of the study was formation of the drops of melt when the heat source was localized at the edge of the wire. Despite the obtained results were approved experimentally in Ref. [8], the approach can be hardy applied to solve the problem in question. However, even if an analytical solution to a heat conduction problem is formally obtained, several problems arise between it and its engineering application. There should be a computational algorithm that transforms the obtained formulae into numerical results, and the algorithm is far from being correct and precise. In order to solve this problem the concept of intrinsic verification of analytical solutions of heat conduction problems found in books or another databases was describes in Ref. [9]. In Ref. [10], the moving mesh method was used to simulate the blowup in a reactionediffusion equation with traveling heat source. It was shown that the finite-time blowup occurred if the speed of the movement of the heat source remained sufficiently low, and the blowup procedure was not fixed at one point not like that for stationary heat source. In this simulation, a new moving mesh
algorithm was designed to deal with the difficulty caused by the delta function in the traveling heat source. The convergence rates were verified and new blowup figures were generated from the numerical experiments. In Ref. [11], an analysis for simulating melting heat transfer around a moving, horizontal, and cylindrical heat source is presented. Motivated by the experimental observations, the melt domain was divided into two regions, namely, the close-contact region and the melt pool region. Two mathematical models were formulated and solution procedures were developed accordingly. The temperature and the flow fields in the two regions are calculated for a constant surface temperature heat source, and the resulting velocity of the source and motion and shape of the interface are determined. The effects of the prescribed surface temperature of the source and its density, as well as influence of natural convection in the melt pool, were investigated and reported. The predicted melt flow structure and the motion and shape of the solideliquid interface are found to be in good agreement with the experimental observations when natural convection in the melt is included in the model. Thus, the melting process in a close proximity of the moving heat source was examined but not in the whole object. The analyses of these and other works show that despite the fact that many solutions for different aspects of the problem were obtained on a high mathematical and physical level, they all are rather far from direct engineering needs. The main tasks of engineering interest are the following. Can a rod be completely melted by a localized stationary heat source? Can this be done by a moving heat source? What program of the heat source motion does allow a complete meltdown of the rod over the shortest period of time? A cell model of heat conduction that formally uses the mathematical tools of the theory of Markov chains to solve the problems is proposed below. It was successfully used to describe the processes in particle technology [12,13] and heat and mass transfer in technological equipment [14,15]. Its application allows freeing the solution from necessity to make fargoing assumptions in order to obtain it that often decreases model adequacy. Its application allows us to get rid of the necessity to make far-going assumptions e which often impairs the model's adequacy. 2. Mathematical formulation of the problem The object of modeling is shown schematically in Fig. 1. It is a one-dimensional rod of the unit cross section, heat insulated at the edges. Its side periphery is fully or partly open to the heat exchange with outside medium. A localized heat source is applied in the point xs than can travel over the rod length. A program of its motion xs(t) is given in advance. If the local rod temperature reaches the melting point tme the phase transformation occurs. It can be melting or solidification depending on the temperature growth or decrease. At
Fig. 1. Schematic presentation of the process.
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
3
the same time the heat emission from the rod to the outside medium occurs. The heat conduction equation for the process can be written in the following form:
cr
vt v vt ¼ l þ qin qout ±qr vt vx vx
(1)
where t(x,t) is the temperature distribution over the rod length at the moment of time t, c is the specific thermal capacity, r is the density, and l is the heat conduction coefficient that is different in different zones of the rod separated by the points x1(t) and x2(t), positions of which are unknown in advance. The initial and boundary conditions for the problem have the form
tðx; 0Þ ¼ t0 ðxÞ
(2) Fig. 2. Schematic presentation of the cell model.
vtð0; tÞ vtðL; tÞ ¼ ¼0 vx vx
(3)
The symbol q is related to the density of heat flow acting on the rod. The localized heat source qin is applied at the point xs(t) that can move over the rod length according to a given in advance program. It can be presented by the formula
qin ðx; tÞ ¼ Qin dðx xs ðtÞÞ
(4)
where Qin is the heat source thermal capacity applied to the rod, d in the symbol of Delta function. If it is supposed that the heat exchange between the rod with the outside medium goes due heat emission, the value of qout can be described as follows
qout ðx; tÞ ¼ aout ðtðx; tÞ tout Þ
Y
(5)
where aout is the heat emission coefficient, tout is the temperature of outside medium, P is the cross-sectional perimeter of the rod. The value qr means the local heat flow density caused by the latent heat of phase transformation. It has the sign minus in the case of melting, and the sign plus in the case of solidification. Thus, in order to describe the process, it is necessary to solve the unsteady non-linear Eq. (1) with the moving boundaries separating the zones with different phase states. Its analytical solution is not realistic particularly because of the non-linearity and moving boundaries inside the rod. As it was mentioned in the Section 1, the authors consider that the appropriate mathematical tool to solve the equation numerically is a cell model of the process that allows not only obtaining an approximate solution but also deriving the basic balance equations of the model on the clear engineering basis. 3. A cell model for numerical solution of the problem The object of modeling is a one-dimensional rod with known thermo-physical properties that is heat insulated at the ends and fully or partly open from the lateral surface. The rod is heated up by a heat source the size of which is much smaller that the size of the rod. The heat source can move over the rod length under this or that program of motion. The zones of the rod that are free from the heat source action can exchange heat with the outside environment of constant temperature. If a local zone of the rod is in the solid state and its temperature rises above the melting point, the process of its melting begins. If the zone is in the liquid state and its temperature falls below the melting point, it leads to its solidification.
The cell model of the process is shown schematically in Fig. 2. The total length of the rod is divided into m cells of the size Dx ¼ L/m, which is small enough to consider all the thermal properties of the cell as being homogeneously distributed over the cell. The thermal state of the rod (i.e., the chain of the cells) can be presented as a set of state column vectors of the size mx1. For example, temperature t and liquid phase content ml distribution can be written as
3 t1 6 t2 7 7 t ¼ tj ¼ 6 4 … 5; tm 2
(6)
2
h
ml ¼ mlj
i
3 ml1 6 ml2 7 7 ¼6 4 … 5 mlm
(7)
where j is the cell's number. The process is observed in discrete moments of time tk ¼ (k 1) Dt where Dt is the time transition duration and k is the time transition number that can be interpreted as the integer-valued analogue of the current time. The vector of distribution of heat over the cells can be calculated as
Q ¼ c: V: r: t
(8)
where c is the vector of specific heat capacities, V is the vector of cell volumes, r is the vector of material densities, and the operator.* means element by element vector multiplication. At each time transition the vector Q changes and this change can be described by the recurrent matrix equation
Q kþ1 ¼ PkQ Q k þ DQ kin þ DQ kout
(9)
where PkQ is the transition matrix of heat conduction that can vary from one transition to another in general case, DQkin is the vector of heat obtained from the heat source during Dt, DQkout is the vector of heat transferred to the outside environment during the same time.
4
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
In order to derive the rule of building the matrix PQ for the rod of variable thermo-physical properties let us calculate the amount of heat that transits from the cell j to the cell j þ 1 during the k-th transition due to the heat conduction
k DQjþ1;j ¼ lkj
¼ Qjk
k Qjþ1
PkQ
lk Dt 6 1 k 1k 6 c1 r1 Dx2 6 6 6 lk1 Dt 6 6 ¼ 6 ck rk Dx2 6 1 1 6 6 6 6 0 4 …
(10)
Dt
where Pkt is the heat conduction matrix for temperature, Dtkin is the vector of temperature variation during one time k transition due to action of a heater, and Dtout is the vector of temperature variation due to the heat exchange with the outside environment. It can be easily shown that the heat conduction matrix for the temperature is, in fact, the transposed matrix for the heat, i.e., Pkt ¼ (PkQ)T. Let us examine the source terms in Eq. (14) and suppose, for the sake of determinacy, that the both heat sources transfer the heat by heat emission. If the heat source is applied to the cell js
0
ck2 rk2 Dx2 lk1 k k c2 r2
Dt Dx
2
lk2
lk2 k k c2 r2
lk2 k k c3 r3
Dt 2
Dx
Dt
ck2 rk2 Dx2 …
1
lk2
Dt Dx2
Dt
2
ck3 rk3 Dx …
The j-th column of the matrix belongs to the cell j. Its element (j þ 1,j) is the part of heat that transits to the next cell due to heat conduction at the k-th transition, the element (j 1,j) is the part of heat that transits to the previous cell, and the element (j,j) is the part of heat that remains in the cell. Each column of the matrix must meet the condition of normalization, i.e., the sum of all elements in each column must be equal to one, and no-one element must be negative. Namely such combination of indices as it shown in Eq. (11) must be used. It concerns the backward transitions from the j-th cell where the value of l must be taken from the previous cell. The point is that the specific heat capacity cj and the density rj are attributes of the cell j but the heat conduction coefficient l is the attribute of a pair of two neighboring cells and must be assigned to but one cell of the two: either to the cell j, or to the cell j þ 1, this choice being arbitrary. The matrix given by Eq. (11) allows describing heat conduction in a composite rod and in a rod with temperature-dependent process parameters, i.e., non-linear heat conduction. In the simplest case of constant thermo-physical properties, the matrix PQ is constant and has the following form
2
PQ
1d 6 d 6 6 0 ¼6 6 … 6 4 0 0
where
d 1 2d d … 0 0
0 d 1 2d … 0 0
(14)
3
lk1
1
(13)
tkþ1 ¼ Pkt tk þ Dtkin þ Dtkout
The first term in the right hand part of the equation can be interpreted as the amount of heat that leaves the cell j due to heat conduction, and the second one as the amount of heat that comes from the cell j þ 1 to the cell j. The expressions in parenthesis are the corresponding parts of the heat accumulated in the cells before the k-th time transition. These parts can be grouped in the tridiagonal transition matrix PQ that has the following form
2
l Dt cr Dx2
The recurrent Eq. (9) written with respect to heat variation can be rewritten with respect to temperature variation
Qjk
SDxrk ck SDxrkjþ1 ckjþ1 j j SDt ¼ lkj SDt Dx Dx ! ! k k lj Dt lj Dt k Qjþ1 2 k k k rj cj Dx rjþ1 ckjþ1 Dx2
k tk tjþ1 j
d¼
… … … … … …
0 0 0 ::: 1 2d d
3 0 0 7 7 0 7 7 ::: 7 7 d 5 1d
(12)
lk3
Dt
ck3 rk3 Dx2
…7 7 7 7 7 7 …7 7 7 7 7 7 …7 5 …
(11)
(the point source in the terms of continuous distribution), all elements of the vector Dtkin are equal to zero except the js-th one, which is
k Dtin ¼ a1 ths tjsk
(15)
where
a1 ¼
a1 S1 Dt : cjs rjs V
(16)
where a1 is the reduced coefficient of heat emission from the source to the rod, a1 is the coefficient of heat emission from the source to the rod, S1 is the surface of heat contact between the source to the rod, V is the cell volume, and ths is the temperature of heat source. The number of the cell js that the heat source is applied to can vary from one transition to another with the residence time Kd (the value of Kd is expressed in the number of transitions, during which the heat source remains at the cell js) according to a given program. k The vector Dtout has the only zero element at j ¼ js. All the other elements can be calculated as
k Dtoutj ¼ a2j tjk tout ; Where
jsjs
(17)
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
a2j ¼
a2 S2 Dt : cj rj V
5
(18)
where a2j is the reduced coefficient of heat emission from the rod to the outside environment, a2j is the coefficient of heat emission from the rod to the outside environment, S2 is the surface of heat contact between the rod and the outside environment, tout is the outside temperature. If the heat source temperature ths is higher than the melting temperature tme, melting of the rod begins in the proximity of the source and propagates along it. In order to describe the phase transformation, it is convenient to introduce the relative mass concentration of the liquid phase Clj ¼ mlj/M where M is the mass of the cell that is supposed to be constant. At each time transition, each cell must be checked for presence or absence of phase transformation (melting or solidification). If the transformation has place, the process must be modeled according to the following algorithm. Melting: if tkj þ 1 > tkj (heating up), tkj þ 1 > tme and Cklj < 1, then k þ 1 tj ¼ tme and
Cljkþ1 ¼ Cljk þ a1 S1 ðths tme ÞDt=r
(19)
where r is the specific latent heat of melting. If Cklj þ 1 > 1 then Cklj þ 1 ¼ 1, i.e., the cell is completely melted and continues to be heated up in the liquid state with different thermal properties in comparison to the solid state. Solidification: if tkj þ 1 < tkj (cooling), tkj þ 1 < tme and Cklj > 0, then k þ 1 tj ¼ tme and
Cljkþ1 ¼ Cljk a1 S1 ðths tme ÞDt=r
(20)
If Cklj þ 1 < 0 then Cklj þ 1 ¼ 0, i.e., the cell is completely solidified and continues to be cooled down in the solid state. Material properties: If Cklj þ 1 ¼ 0 then lkj þ 1 ¼ ls, ckj þ 1 ¼ cs; if Cklj þ 1 ¼ 1 then lkj þ 1 ¼ ll, ckj þ 1 ¼ cl. Thus, the cell model presented above allows describing the nonlinear heat conduction in a rod accompanied by the phase transformation melting/solidification. As far as the cells models are a version of the finite difference method, its accuracy can be estimated on the same basis as the method. If Dx and Dt limit to zero the numerical solution obtained with a cell model limits to the exact solution to Eq. (1). The choice of Dx and Dt can be done in numerical experiments when the solution becomes Dx and Dt independent with a required accuracy.
4. Results and discussion This section presents some of the results of numerical experiments with the model and their discussion. The following parameters for modeling were used: ls ¼ 36 W/(m·K), cs ¼ 130 J/ (kg·K), ll ¼ 34 W/(m·K), cl ¼ 150 J/(kg·K), rs ¼ rl ¼ 104 kg/m3, r ¼ 25·103 J/kg, a1 ¼ 0.5, a2 ¼ 0.015. All the temperatures are presented in the dimensionless form (divided by the heat source temperature) so that ths ¼ 1, tme ¼ 0.5, tout ¼ 0. The initial temperature of the rod is taken zero. The number of cells is taken as being equal to 11 in order to make pictorial presentation of phase transformation more obvious and visually demonstrate the main features that appear at different programs of the heat source motion. In order to explain the physical features of the process expressed by the cell model let us examine the evolution of the thermal state of the rod when the heat source is immovable and
Fig. 3. Evolution of temperature distribution (a) and relative content of melt (b) for immovable heat source.
applied to the middle of the rod, i.e., to the cell 6. The results are shown in Fig. 3. First, the central cell is heated up to the melting point. At the same time other cells are heated from the central one due to heat conduction. When the temperature of the central cell reaches the melting point melting of it begins at the constant temperature tme. Melting of neighboring cells cannot begin until melting of the central cell is completed. Only after that the cell can increase its temperature due to heating it in the liquid state and heat up the neighboring cells to the melting point and so on. It is more visual to present the results of solid/liquid state evolution as contour graphs that is done in the next figures. Fig. 4 shows the melting front propagation at different conditions of heating the rod. The upper graph is related to the fixedposition heat source localized in the central cell of the rod that was examined above. The waviness of the lines is conditioned by the fact that if the phase transformation occurs in a cell and its temperature is constant, the next cell cannot be heated up to the melting temperature due to heat conduction and melting does not begin in the cell until the neighboring cell is completely melted. If the rod is divided into a larger number of cells, the zone of the solid/ liquid mixture becomes narrower but numerical experiments show that the centerline remains approximately the same. It can be seen from the graph that the rod cannot be melted completely: the very edges of it are only melted partly. The graph in the middle shows the propagation of the melting front for the two-point program of the heat source motion js ¼ 4, 8, 4, … at the residence time of the source at each cell of the two equal
6
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
Fig. 4. Propagation of the melting front at various residence time of the heat source at the cells for two-point program of the source motion.
to 2 time transitions (Kd ¼ 2). In this case the melting front moves much faster, and after 280 time transitions the rod appears to be completely melted. Switching to the residence time equal to 4 slows the process down to a certain extent but not considerably. It is necessary to note that at a really long residence time the complete meltdown of the rod becomes problematic again because the solidified part of the rod appears periodically at one or the other end of the rod. Fig. 5 shows how the source positions at the two-point program of its motion influence the melting front propagation. If the points of the source application are very close to the rod's center the process of melting the rod completely takes a much longer time (upper graph) in comparison to the case when they are close to the middles of the rod halves (middle graph). The longest time is required when the rod is heated at its very edges (lower graph). These results are generalized in the graph shown in Fig. 6 where the objective function is the time (number of time transitions) Kcm required for the complete meltdown of the rod. The positions of the heat source application are presented as the distance Z from the rod center. The value Kd ¼ 1 is not included in the graph because, under this residence time, complete meltdown of the rod is impossible at any program of the heat source motion. It can be seen from the graph that the rod can be completely melted during the shortest time at Kd ¼ 2 and Z ¼ 3, i.e., at js ¼ 3, 9, 3,…, the growth of Kd has practically no influence at that.
Fig. 5. Propagation of the melting front at various positions of the heat source for twopoints program of the source motion (Kd ¼ 2).
It is necessary to note that the graph shown in Fig. 6 is not the universal result of optimization and strongly depends on material properties and dimensions of the rod. For instance, if we take a rod of the same kind but twice as long, there is no program of the heat source motion that allows for its complete meltdown. In order to do that, we would have either to increase the heat source temperature, or to decrease heat emission to the outside environment. The latter is hardly realistic.
Fig. 6. Duration of complete melting the rod versus the residence time of the heat source in the cells and their distance from the rod center.
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
7
Acknowledgments This work was supported by the Russian Foundation for Basic Research (project 15-08-01684). Nomenclature a1, a2
Fig. 7. Propagation of the melting front at decreasing temperature of the heat source for the two-points program of the source motion (Z ¼ 2, Kd ¼ 4).
In order to demonstrate the possibilities of the model, the graph of melting front propagation at the source temperature decreasing with time for the two-point program is shown in Fig. 7. At the first stage of the process melting goes similar to one shown in the middle graph of Fig. 4, then melting is getting stabilized, and then the inverse solidification begins that ends after 750 time transitions. It can be easily seen that the proposed model and algorithm of calculations can be applied to even more complex cases of heating. For instance, the convection heating can be replaced by induction heating e it is enough to change Eqs. 15 and 16 with the equations for induction heat transfer. The heat source localized at a single cell can be replaced by a moving heat source distributed over several cells. And finally, the model can be generalized to melting a rectangular plate, or a plate of arbitrary shape. However, despite the fact that the plate model is based on the same principles, special work should be done to build the transition matrix and other model variables, which is beyond the objective of the present study. 5. Conclusions A cell model to describe melting a rod by a localized heat source travelling over the rod length due to specially chosen program is proposed. The model allows predicting the evolution of temperature distribution and melting front propagation for different programs of heat source motion. It is shown, that when melting the rod completely by a fixed-position heat source is impossible, moving the heat source in accordance with a special two-point program allows for achieving its complete meltdown. The melting rate strongly depends on the positions of the cells of heat source application and on the residence time of the source at the cells of application. The optimal program of melting the rod with taken properties is found that allows for its complete meltdown that occurs during the shortest time. It is demonstrated how the model works under the variable environment temperature. Some of the directions of the further development of this model are discussed.
Cl cj, c j k Kd Kcm m ml, ml PQ Pt Qj, Q q r S1, S2 tj, t Dx Z
dimensionless parameters of heat transfer to and from the rod (e) relative mass concentration of liquid phase (e) specific heat capacity (J kg1 K1) cell number (e) transition number (e) residence time (in the number of transitions) (e) number of transitions for a complete meltdown (e) number of cells (e) mass of liquid phase (kg) matrix of heat conduction for heat (e) matrix of heat conduction for temperature (e) heat (J) Heat flow density (Jm1s1) specific latent heat of melting/solidification (Jkg1) area of heat transfer to and from the rod (m2) temperature ( C) cell length (m) distance from the central cell of the rod (in cells number) (e)
Greek symbols a heat transfer coefficient (Wm2 s1 K1) Dt transition duration (s) l heat conduction coefficient (Wm1K1) r density (kg m3) Subscripts hs heat source l liquid phase me melting s solid phase out outside environment References [1] J. Kidawa-Kukla, Temperature distribution in a rectangular plate heated by a moving heat source, Int. J. Heat Mass Transf. 51 (2008) 865e872. [2] M. Shibahara, S.N. Atluri, The meshless local Petrov-Galerkin method for the analysis of heat conduction due to a moving heat source, in welding, Int. J. Therm. Sci. 50 (2011) 984e992. [3] M. Van Elsen, M. Baelmans, P. Mercelis, J.-P. Kruth, Solutions for modeling moving heat source in a semi-infinite medium and applications to laser material processing, Int. J. Heat Mass Transf. 50 (2007) 4872e4882. [4] H. Belghazi, M. El Ganaoui, J.C. Labbe, Analytical solution of unsteady heat conduction in a two-layered material in imperfect contact subjected to a moving heat source, Int. J. Therm. Sci. 49 (2010) 311e318. [5] J. Winczek, Analytical solution to transient temperature field in a half-infinite body caused by moving volumetric heat source, Int. J. Heat Mass Transf. 53 (2010) 5774e5781. [6] J. Winczek, New approach to modeling of temperature field in surfaced steel elements, Int. J. Heat Mass Transf. 54 (2011) 4702e4709. [7] L.J. Huang, P.S. Ayyaswamy, I.M. Cohen, Melting and solidification of thin wires: a class of phase-change problems with a mobile interface e I. Analysis, Int. J. Heat Mass Transf. 38 (1995) 1637e1645. [8] I.M. Cohen, L.J. Huang, P.S. Ayyaswamy, Melting and solidification o thin wires: class of phase-change problems with a mobile interface e II. Experimental confirmation, Int. J. Heat Mass Transf. 38 (1995) 1647e1659. [9] K.D. Cole, J.V. Beck, K.A. Woodbury, F. de Monte, Intrinsic verification and a heat conduction database, Int. J. Therm. Sci. 78 (2014) 36e47. [10] Jingtang Ma, Yingjun Jiang, Moving mesh methods for blowup in reactionediffusion equations with traveling heat source, J. Comput. Phys. 228 (2009) 6977e6990.
8
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
[11] M.K. Moallemi, R. Visksn, Analysis of melting around a moving heat source, Int. J. Heat Mass Trans. 29 (1986) 1271e1282. [12] H. Berthiaux, V. Mizonov, Applications of Markov chains in particulate process engineering: a review, Can. J. Chem. Eng. 85 (2004) 1143e1168. [13] H. Berthiaux, V. Mizonov, V. Zhukov, Application of the theory of Markov chains to model different processes in particle technology, Powder Technol. 157 (2005) 128e137.
[14] V. Mizonov, H. Berthiaux, P. Arlabosse, D. Djerroud, Application of the theory of Markov chains to model heat and mass transfer between stochastically moving particulate and gas flows, Granul. Matter 10 (2008) 335e340. [15] V. Mizonov, N. Yelin, P. Yakimychev, A cell model to describe and optimize heat and mass transfer in contact heat exchangers, Energy Power Eng. 3 (2011) 144e149.
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Numerical study of melting a rod by a periodically moving local heat source Vadim Mizonov a, *, Nicolay Yelin b, 1 a b
Department of Applied Mathematics, Ivanovo State Power Engineering University, Rabfakovskaya 34, 153003, Ivanovo, Russia Department of Hydraulics, Heat and Mass Transfers, Ivanovo State Polytechnic University, 8 Marta 20, 153037, Ivanovo, Russia
a r t i c l e i n f o
a b s t r a c t
Article history: Received 21 July 2014 Received in revised form 14 June 2015 Accepted 15 June 2015 Available online xxx
The objective of the study is developing a simple yet informative mathematical model that describes the kinetics of melting a rod under the action of a localized periodically moving heat source. For this purpose, a cell model is used with the heat conduction matrix that takes into account different properties of liquid and solid phases of the rod material. Zones of the rod that are outside of the local heat source action have heat exchange with the outside environment. It is shown that the melting kinetics strongly depends on the program of heat source motion along the rod and on its residence time at each of its positions. The optimal program of heat source motion that allows melting the whole rod over the shortest possible period of time is found. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Heat conduction Melting Solidification Local heat source Heat transfer Optimal program
1. Introduction The problem in question arises in some technologies when there is a need to heat up and melt an object by a heat source whose dimensions are much smaller than those of the object to be melted. The simplest example is defrosting a finite fragment of a water pipe full of frozen water by a burner. Sometimes (for instance, at a low outside temperature), this cannot be done if the heat source is only applied to one particular point of the fragment, and it is only its motion along the pipe that can solve the problem. In this case, the complete defrosting time strongly depends on the program of the burned motion over the pipe length e which can be the objective function of investigation. Analogous problems can arise in other technologies, for example, in hand treatment of roof materials by a gas burner, in soldering, and in additive technologies when laser beam treatment of thin layer of fine powder occurs. This problem is related to the non-linear heat conduction problems with a localized moving heat source. Its non-linearity is
* Corresponding author. Tel.: þ7 910 9948858; fax: þ7 4932385701. E-mail addresses: [email protected] (V. Mizonov), [email protected] (N. Yelin). 1 Tel.: þ7 910 9810020. http://dx.doi.org/10.1016/j.ijthermalsci.2015.06.005 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.
conditioned by the phase transformation of melting/solidification type and various thermal effects connected with it. Different aspects of the problem solution can be found in literature. However, most of them are connected with the analytical or semi-analytical approach to solving this problem e which inevitably implies using far-going assumptions. A solution to the problem of heat conduction in a rectangular plate exposed to a moving heat source was presented in Ref. [1]. The heat source moved along an elliptical trajectory that always remained within the boundaries of the plate area. The exact solution to the problem in an analytical form was obtained by applying the Green's function method. Exemplary results of numerical calculations to determine the temperature distribution in the plate were presented. However, first of all, this was the only trajectory of the source motion and, second, the phase transformation was not taken into account. The problem of the object being heated by a moving source with application to welding was also examined in Ref. [2] where the meshless local Petrov-Galerkin method was developed. The problem of bead-on-plate welding (a moving heat source problem), and the fundamental properties of the method were investigated to verify the applicability of the proposed method. The investigation was mainly oriented to the influence of the number of nodal points on the accuracy of the solution. The issue of the influence of the program of the source motion on the
2
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
process properties was beyond this particular study. The study presented in Ref. [3] describes the analytical and numerical solution of the heat conduction equation for a localized moving heat source of any type to be used in laser material processing, such as welding, layered manufacturing and laser alloying. The analytical solution for a uniform heat source was derived from the solution of an instantaneous point heat source. The result was evaluated numerically and was compared to existing solutions for the moving point source and a semi-ellipsoidal source. Next, the result was used to demonstrate how such model can be used to study the effect of the heat source geometry. To investigate the effect of the temperature dependence of the material parameters, in particular the latent heat of fusion, a finite difference model was implemented. In spite of the fact that many components of the process were taken into account the influence of the program of the source motion remained beyond the scope of attention of the study. An analytical approach of transient heat conduction in two-layered material, of finite depth, with an imperfect thermal contact, subjected to a moving gaussian laser beam was developed in Ref. [4]. The method consists in deriving the solution of the homogeneous part of the heat equation by using the well known separation of variables method and expressing the source term in series form. This model has been successfully applied on a practical system; laser cladding of electronic copper tracks on alumina substrates. This analytical model can also be used for estimation of the thermal contact resistance between the layers. However, the model does not take into account the phase transformation. An analytical method of computation of temperature field in a half-infinite body caused by heat source with changeable direction of motion was presented in Ref. [5]. Analytical temperature field was approximated by straight segments for volumetric heat source with a trajectory considering temperature changes caused by next transitions (increase in temperature connected to action of the heat source and self-cooling of areas heated-up earlier). In that instance computations were carried out for cuboidal elements made of steel for various heat trajectories. However, the phase transformation in the body was not taken into account too. Later on, in the paper [6], the melting process was introduced in the model but only as a reason of the increase of heat transfer from the heat source to the body. Melting of the body itself was not examined. A model that describes the transient heating of a thin wire causing the tip to melt, roll-up of the molten mass into a ball due to surface tension forces, and the subsequent solidification of the molten material due to conduction up the wire and convection and radiation from the surface, was proposed in Ref. [7]. The wire was assumed to be heated at its lower tip to a temperature beyond the melting temperature of the wire material by heat flux from an electrical discharge. However, the objective of the study was formation of the drops of melt when the heat source was localized at the edge of the wire. Despite the obtained results were approved experimentally in Ref. [8], the approach can be hardy applied to solve the problem in question. However, even if an analytical solution to a heat conduction problem is formally obtained, several problems arise between it and its engineering application. There should be a computational algorithm that transforms the obtained formulae into numerical results, and the algorithm is far from being correct and precise. In order to solve this problem the concept of intrinsic verification of analytical solutions of heat conduction problems found in books or another databases was describes in Ref. [9]. In Ref. [10], the moving mesh method was used to simulate the blowup in a reactionediffusion equation with traveling heat source. It was shown that the finite-time blowup occurred if the speed of the movement of the heat source remained sufficiently low, and the blowup procedure was not fixed at one point not like that for stationary heat source. In this simulation, a new moving mesh
algorithm was designed to deal with the difficulty caused by the delta function in the traveling heat source. The convergence rates were verified and new blowup figures were generated from the numerical experiments. In Ref. [11], an analysis for simulating melting heat transfer around a moving, horizontal, and cylindrical heat source is presented. Motivated by the experimental observations, the melt domain was divided into two regions, namely, the close-contact region and the melt pool region. Two mathematical models were formulated and solution procedures were developed accordingly. The temperature and the flow fields in the two regions are calculated for a constant surface temperature heat source, and the resulting velocity of the source and motion and shape of the interface are determined. The effects of the prescribed surface temperature of the source and its density, as well as influence of natural convection in the melt pool, were investigated and reported. The predicted melt flow structure and the motion and shape of the solideliquid interface are found to be in good agreement with the experimental observations when natural convection in the melt is included in the model. Thus, the melting process in a close proximity of the moving heat source was examined but not in the whole object. The analyses of these and other works show that despite the fact that many solutions for different aspects of the problem were obtained on a high mathematical and physical level, they all are rather far from direct engineering needs. The main tasks of engineering interest are the following. Can a rod be completely melted by a localized stationary heat source? Can this be done by a moving heat source? What program of the heat source motion does allow a complete meltdown of the rod over the shortest period of time? A cell model of heat conduction that formally uses the mathematical tools of the theory of Markov chains to solve the problems is proposed below. It was successfully used to describe the processes in particle technology [12,13] and heat and mass transfer in technological equipment [14,15]. Its application allows freeing the solution from necessity to make fargoing assumptions in order to obtain it that often decreases model adequacy. Its application allows us to get rid of the necessity to make far-going assumptions e which often impairs the model's adequacy. 2. Mathematical formulation of the problem The object of modeling is shown schematically in Fig. 1. It is a one-dimensional rod of the unit cross section, heat insulated at the edges. Its side periphery is fully or partly open to the heat exchange with outside medium. A localized heat source is applied in the point xs than can travel over the rod length. A program of its motion xs(t) is given in advance. If the local rod temperature reaches the melting point tme the phase transformation occurs. It can be melting or solidification depending on the temperature growth or decrease. At
Fig. 1. Schematic presentation of the process.
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
3
the same time the heat emission from the rod to the outside medium occurs. The heat conduction equation for the process can be written in the following form:
cr
vt v vt ¼ l þ qin qout ±qr vt vx vx
(1)
where t(x,t) is the temperature distribution over the rod length at the moment of time t, c is the specific thermal capacity, r is the density, and l is the heat conduction coefficient that is different in different zones of the rod separated by the points x1(t) and x2(t), positions of which are unknown in advance. The initial and boundary conditions for the problem have the form
tðx; 0Þ ¼ t0 ðxÞ
(2) Fig. 2. Schematic presentation of the cell model.
vtð0; tÞ vtðL; tÞ ¼ ¼0 vx vx
(3)
The symbol q is related to the density of heat flow acting on the rod. The localized heat source qin is applied at the point xs(t) that can move over the rod length according to a given in advance program. It can be presented by the formula
qin ðx; tÞ ¼ Qin dðx xs ðtÞÞ
(4)
where Qin is the heat source thermal capacity applied to the rod, d in the symbol of Delta function. If it is supposed that the heat exchange between the rod with the outside medium goes due heat emission, the value of qout can be described as follows
qout ðx; tÞ ¼ aout ðtðx; tÞ tout Þ
Y
(5)
where aout is the heat emission coefficient, tout is the temperature of outside medium, P is the cross-sectional perimeter of the rod. The value qr means the local heat flow density caused by the latent heat of phase transformation. It has the sign minus in the case of melting, and the sign plus in the case of solidification. Thus, in order to describe the process, it is necessary to solve the unsteady non-linear Eq. (1) with the moving boundaries separating the zones with different phase states. Its analytical solution is not realistic particularly because of the non-linearity and moving boundaries inside the rod. As it was mentioned in the Section 1, the authors consider that the appropriate mathematical tool to solve the equation numerically is a cell model of the process that allows not only obtaining an approximate solution but also deriving the basic balance equations of the model on the clear engineering basis. 3. A cell model for numerical solution of the problem The object of modeling is a one-dimensional rod with known thermo-physical properties that is heat insulated at the ends and fully or partly open from the lateral surface. The rod is heated up by a heat source the size of which is much smaller that the size of the rod. The heat source can move over the rod length under this or that program of motion. The zones of the rod that are free from the heat source action can exchange heat with the outside environment of constant temperature. If a local zone of the rod is in the solid state and its temperature rises above the melting point, the process of its melting begins. If the zone is in the liquid state and its temperature falls below the melting point, it leads to its solidification.
The cell model of the process is shown schematically in Fig. 2. The total length of the rod is divided into m cells of the size Dx ¼ L/m, which is small enough to consider all the thermal properties of the cell as being homogeneously distributed over the cell. The thermal state of the rod (i.e., the chain of the cells) can be presented as a set of state column vectors of the size mx1. For example, temperature t and liquid phase content ml distribution can be written as
3 t1 6 t2 7 7 t ¼ tj ¼ 6 4 … 5; tm 2
(6)
2
h
ml ¼ mlj
i
3 ml1 6 ml2 7 7 ¼6 4 … 5 mlm
(7)
where j is the cell's number. The process is observed in discrete moments of time tk ¼ (k 1) Dt where Dt is the time transition duration and k is the time transition number that can be interpreted as the integer-valued analogue of the current time. The vector of distribution of heat over the cells can be calculated as
Q ¼ c: V: r: t
(8)
where c is the vector of specific heat capacities, V is the vector of cell volumes, r is the vector of material densities, and the operator.* means element by element vector multiplication. At each time transition the vector Q changes and this change can be described by the recurrent matrix equation
Q kþ1 ¼ PkQ Q k þ DQ kin þ DQ kout
(9)
where PkQ is the transition matrix of heat conduction that can vary from one transition to another in general case, DQkin is the vector of heat obtained from the heat source during Dt, DQkout is the vector of heat transferred to the outside environment during the same time.
4
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
In order to derive the rule of building the matrix PQ for the rod of variable thermo-physical properties let us calculate the amount of heat that transits from the cell j to the cell j þ 1 during the k-th transition due to the heat conduction
k DQjþ1;j ¼ lkj
¼ Qjk
k Qjþ1
PkQ
lk Dt 6 1 k 1k 6 c1 r1 Dx2 6 6 6 lk1 Dt 6 6 ¼ 6 ck rk Dx2 6 1 1 6 6 6 6 0 4 …
(10)
Dt
where Pkt is the heat conduction matrix for temperature, Dtkin is the vector of temperature variation during one time k transition due to action of a heater, and Dtout is the vector of temperature variation due to the heat exchange with the outside environment. It can be easily shown that the heat conduction matrix for the temperature is, in fact, the transposed matrix for the heat, i.e., Pkt ¼ (PkQ)T. Let us examine the source terms in Eq. (14) and suppose, for the sake of determinacy, that the both heat sources transfer the heat by heat emission. If the heat source is applied to the cell js
0
ck2 rk2 Dx2 lk1 k k c2 r2
Dt Dx
2
lk2
lk2 k k c2 r2
lk2 k k c3 r3
Dt 2
Dx
Dt
ck2 rk2 Dx2 …
1
lk2
Dt Dx2
Dt
2
ck3 rk3 Dx …
The j-th column of the matrix belongs to the cell j. Its element (j þ 1,j) is the part of heat that transits to the next cell due to heat conduction at the k-th transition, the element (j 1,j) is the part of heat that transits to the previous cell, and the element (j,j) is the part of heat that remains in the cell. Each column of the matrix must meet the condition of normalization, i.e., the sum of all elements in each column must be equal to one, and no-one element must be negative. Namely such combination of indices as it shown in Eq. (11) must be used. It concerns the backward transitions from the j-th cell where the value of l must be taken from the previous cell. The point is that the specific heat capacity cj and the density rj are attributes of the cell j but the heat conduction coefficient l is the attribute of a pair of two neighboring cells and must be assigned to but one cell of the two: either to the cell j, or to the cell j þ 1, this choice being arbitrary. The matrix given by Eq. (11) allows describing heat conduction in a composite rod and in a rod with temperature-dependent process parameters, i.e., non-linear heat conduction. In the simplest case of constant thermo-physical properties, the matrix PQ is constant and has the following form
2
PQ
1d 6 d 6 6 0 ¼6 6 … 6 4 0 0
where
d 1 2d d … 0 0
0 d 1 2d … 0 0
(14)
3
lk1
1
(13)
tkþ1 ¼ Pkt tk þ Dtkin þ Dtkout
The first term in the right hand part of the equation can be interpreted as the amount of heat that leaves the cell j due to heat conduction, and the second one as the amount of heat that comes from the cell j þ 1 to the cell j. The expressions in parenthesis are the corresponding parts of the heat accumulated in the cells before the k-th time transition. These parts can be grouped in the tridiagonal transition matrix PQ that has the following form
2
l Dt cr Dx2
The recurrent Eq. (9) written with respect to heat variation can be rewritten with respect to temperature variation
Qjk
SDxrk ck SDxrkjþ1 ckjþ1 j j SDt ¼ lkj SDt Dx Dx ! ! k k lj Dt lj Dt k Qjþ1 2 k k k rj cj Dx rjþ1 ckjþ1 Dx2
k tk tjþ1 j
d¼
… … … … … …
0 0 0 ::: 1 2d d
3 0 0 7 7 0 7 7 ::: 7 7 d 5 1d
(12)
lk3
Dt
ck3 rk3 Dx2
…7 7 7 7 7 7 …7 7 7 7 7 7 …7 5 …
(11)
(the point source in the terms of continuous distribution), all elements of the vector Dtkin are equal to zero except the js-th one, which is
k Dtin ¼ a1 ths tjsk
(15)
where
a1 ¼
a1 S1 Dt : cjs rjs V
(16)
where a1 is the reduced coefficient of heat emission from the source to the rod, a1 is the coefficient of heat emission from the source to the rod, S1 is the surface of heat contact between the source to the rod, V is the cell volume, and ths is the temperature of heat source. The number of the cell js that the heat source is applied to can vary from one transition to another with the residence time Kd (the value of Kd is expressed in the number of transitions, during which the heat source remains at the cell js) according to a given program. k The vector Dtout has the only zero element at j ¼ js. All the other elements can be calculated as
k Dtoutj ¼ a2j tjk tout ; Where
jsjs
(17)
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
a2j ¼
a2 S2 Dt : cj rj V
5
(18)
where a2j is the reduced coefficient of heat emission from the rod to the outside environment, a2j is the coefficient of heat emission from the rod to the outside environment, S2 is the surface of heat contact between the rod and the outside environment, tout is the outside temperature. If the heat source temperature ths is higher than the melting temperature tme, melting of the rod begins in the proximity of the source and propagates along it. In order to describe the phase transformation, it is convenient to introduce the relative mass concentration of the liquid phase Clj ¼ mlj/M where M is the mass of the cell that is supposed to be constant. At each time transition, each cell must be checked for presence or absence of phase transformation (melting or solidification). If the transformation has place, the process must be modeled according to the following algorithm. Melting: if tkj þ 1 > tkj (heating up), tkj þ 1 > tme and Cklj < 1, then k þ 1 tj ¼ tme and
Cljkþ1 ¼ Cljk þ a1 S1 ðths tme ÞDt=r
(19)
where r is the specific latent heat of melting. If Cklj þ 1 > 1 then Cklj þ 1 ¼ 1, i.e., the cell is completely melted and continues to be heated up in the liquid state with different thermal properties in comparison to the solid state. Solidification: if tkj þ 1 < tkj (cooling), tkj þ 1 < tme and Cklj > 0, then k þ 1 tj ¼ tme and
Cljkþ1 ¼ Cljk a1 S1 ðths tme ÞDt=r
(20)
If Cklj þ 1 < 0 then Cklj þ 1 ¼ 0, i.e., the cell is completely solidified and continues to be cooled down in the solid state. Material properties: If Cklj þ 1 ¼ 0 then lkj þ 1 ¼ ls, ckj þ 1 ¼ cs; if Cklj þ 1 ¼ 1 then lkj þ 1 ¼ ll, ckj þ 1 ¼ cl. Thus, the cell model presented above allows describing the nonlinear heat conduction in a rod accompanied by the phase transformation melting/solidification. As far as the cells models are a version of the finite difference method, its accuracy can be estimated on the same basis as the method. If Dx and Dt limit to zero the numerical solution obtained with a cell model limits to the exact solution to Eq. (1). The choice of Dx and Dt can be done in numerical experiments when the solution becomes Dx and Dt independent with a required accuracy.
4. Results and discussion This section presents some of the results of numerical experiments with the model and their discussion. The following parameters for modeling were used: ls ¼ 36 W/(m·K), cs ¼ 130 J/ (kg·K), ll ¼ 34 W/(m·K), cl ¼ 150 J/(kg·K), rs ¼ rl ¼ 104 kg/m3, r ¼ 25·103 J/kg, a1 ¼ 0.5, a2 ¼ 0.015. All the temperatures are presented in the dimensionless form (divided by the heat source temperature) so that ths ¼ 1, tme ¼ 0.5, tout ¼ 0. The initial temperature of the rod is taken zero. The number of cells is taken as being equal to 11 in order to make pictorial presentation of phase transformation more obvious and visually demonstrate the main features that appear at different programs of the heat source motion. In order to explain the physical features of the process expressed by the cell model let us examine the evolution of the thermal state of the rod when the heat source is immovable and
Fig. 3. Evolution of temperature distribution (a) and relative content of melt (b) for immovable heat source.
applied to the middle of the rod, i.e., to the cell 6. The results are shown in Fig. 3. First, the central cell is heated up to the melting point. At the same time other cells are heated from the central one due to heat conduction. When the temperature of the central cell reaches the melting point melting of it begins at the constant temperature tme. Melting of neighboring cells cannot begin until melting of the central cell is completed. Only after that the cell can increase its temperature due to heating it in the liquid state and heat up the neighboring cells to the melting point and so on. It is more visual to present the results of solid/liquid state evolution as contour graphs that is done in the next figures. Fig. 4 shows the melting front propagation at different conditions of heating the rod. The upper graph is related to the fixedposition heat source localized in the central cell of the rod that was examined above. The waviness of the lines is conditioned by the fact that if the phase transformation occurs in a cell and its temperature is constant, the next cell cannot be heated up to the melting temperature due to heat conduction and melting does not begin in the cell until the neighboring cell is completely melted. If the rod is divided into a larger number of cells, the zone of the solid/ liquid mixture becomes narrower but numerical experiments show that the centerline remains approximately the same. It can be seen from the graph that the rod cannot be melted completely: the very edges of it are only melted partly. The graph in the middle shows the propagation of the melting front for the two-point program of the heat source motion js ¼ 4, 8, 4, … at the residence time of the source at each cell of the two equal
6
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
Fig. 4. Propagation of the melting front at various residence time of the heat source at the cells for two-point program of the source motion.
to 2 time transitions (Kd ¼ 2). In this case the melting front moves much faster, and after 280 time transitions the rod appears to be completely melted. Switching to the residence time equal to 4 slows the process down to a certain extent but not considerably. It is necessary to note that at a really long residence time the complete meltdown of the rod becomes problematic again because the solidified part of the rod appears periodically at one or the other end of the rod. Fig. 5 shows how the source positions at the two-point program of its motion influence the melting front propagation. If the points of the source application are very close to the rod's center the process of melting the rod completely takes a much longer time (upper graph) in comparison to the case when they are close to the middles of the rod halves (middle graph). The longest time is required when the rod is heated at its very edges (lower graph). These results are generalized in the graph shown in Fig. 6 where the objective function is the time (number of time transitions) Kcm required for the complete meltdown of the rod. The positions of the heat source application are presented as the distance Z from the rod center. The value Kd ¼ 1 is not included in the graph because, under this residence time, complete meltdown of the rod is impossible at any program of the heat source motion. It can be seen from the graph that the rod can be completely melted during the shortest time at Kd ¼ 2 and Z ¼ 3, i.e., at js ¼ 3, 9, 3,…, the growth of Kd has practically no influence at that.
Fig. 5. Propagation of the melting front at various positions of the heat source for twopoints program of the source motion (Kd ¼ 2).
It is necessary to note that the graph shown in Fig. 6 is not the universal result of optimization and strongly depends on material properties and dimensions of the rod. For instance, if we take a rod of the same kind but twice as long, there is no program of the heat source motion that allows for its complete meltdown. In order to do that, we would have either to increase the heat source temperature, or to decrease heat emission to the outside environment. The latter is hardly realistic.
Fig. 6. Duration of complete melting the rod versus the residence time of the heat source in the cells and their distance from the rod center.
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
7
Acknowledgments This work was supported by the Russian Foundation for Basic Research (project 15-08-01684). Nomenclature a1, a2
Fig. 7. Propagation of the melting front at decreasing temperature of the heat source for the two-points program of the source motion (Z ¼ 2, Kd ¼ 4).
In order to demonstrate the possibilities of the model, the graph of melting front propagation at the source temperature decreasing with time for the two-point program is shown in Fig. 7. At the first stage of the process melting goes similar to one shown in the middle graph of Fig. 4, then melting is getting stabilized, and then the inverse solidification begins that ends after 750 time transitions. It can be easily seen that the proposed model and algorithm of calculations can be applied to even more complex cases of heating. For instance, the convection heating can be replaced by induction heating e it is enough to change Eqs. 15 and 16 with the equations for induction heat transfer. The heat source localized at a single cell can be replaced by a moving heat source distributed over several cells. And finally, the model can be generalized to melting a rectangular plate, or a plate of arbitrary shape. However, despite the fact that the plate model is based on the same principles, special work should be done to build the transition matrix and other model variables, which is beyond the objective of the present study. 5. Conclusions A cell model to describe melting a rod by a localized heat source travelling over the rod length due to specially chosen program is proposed. The model allows predicting the evolution of temperature distribution and melting front propagation for different programs of heat source motion. It is shown, that when melting the rod completely by a fixed-position heat source is impossible, moving the heat source in accordance with a special two-point program allows for achieving its complete meltdown. The melting rate strongly depends on the positions of the cells of heat source application and on the residence time of the source at the cells of application. The optimal program of melting the rod with taken properties is found that allows for its complete meltdown that occurs during the shortest time. It is demonstrated how the model works under the variable environment temperature. Some of the directions of the further development of this model are discussed.
Cl cj, c j k Kd Kcm m ml, ml PQ Pt Qj, Q q r S1, S2 tj, t Dx Z
dimensionless parameters of heat transfer to and from the rod (e) relative mass concentration of liquid phase (e) specific heat capacity (J kg1 K1) cell number (e) transition number (e) residence time (in the number of transitions) (e) number of transitions for a complete meltdown (e) number of cells (e) mass of liquid phase (kg) matrix of heat conduction for heat (e) matrix of heat conduction for temperature (e) heat (J) Heat flow density (Jm1s1) specific latent heat of melting/solidification (Jkg1) area of heat transfer to and from the rod (m2) temperature ( C) cell length (m) distance from the central cell of the rod (in cells number) (e)
Greek symbols a heat transfer coefficient (Wm2 s1 K1) Dt transition duration (s) l heat conduction coefficient (Wm1K1) r density (kg m3) Subscripts hs heat source l liquid phase me melting s solid phase out outside environment References [1] J. Kidawa-Kukla, Temperature distribution in a rectangular plate heated by a moving heat source, Int. J. Heat Mass Transf. 51 (2008) 865e872. [2] M. Shibahara, S.N. Atluri, The meshless local Petrov-Galerkin method for the analysis of heat conduction due to a moving heat source, in welding, Int. J. Therm. Sci. 50 (2011) 984e992. [3] M. Van Elsen, M. Baelmans, P. Mercelis, J.-P. Kruth, Solutions for modeling moving heat source in a semi-infinite medium and applications to laser material processing, Int. J. Heat Mass Transf. 50 (2007) 4872e4882. [4] H. Belghazi, M. El Ganaoui, J.C. Labbe, Analytical solution of unsteady heat conduction in a two-layered material in imperfect contact subjected to a moving heat source, Int. J. Therm. Sci. 49 (2010) 311e318. [5] J. Winczek, Analytical solution to transient temperature field in a half-infinite body caused by moving volumetric heat source, Int. J. Heat Mass Transf. 53 (2010) 5774e5781. [6] J. Winczek, New approach to modeling of temperature field in surfaced steel elements, Int. J. Heat Mass Transf. 54 (2011) 4702e4709. [7] L.J. Huang, P.S. Ayyaswamy, I.M. Cohen, Melting and solidification of thin wires: a class of phase-change problems with a mobile interface e I. Analysis, Int. J. Heat Mass Transf. 38 (1995) 1637e1645. [8] I.M. Cohen, L.J. Huang, P.S. Ayyaswamy, Melting and solidification o thin wires: class of phase-change problems with a mobile interface e II. Experimental confirmation, Int. J. Heat Mass Transf. 38 (1995) 1647e1659. [9] K.D. Cole, J.V. Beck, K.A. Woodbury, F. de Monte, Intrinsic verification and a heat conduction database, Int. J. Therm. Sci. 78 (2014) 36e47. [10] Jingtang Ma, Yingjun Jiang, Moving mesh methods for blowup in reactionediffusion equations with traveling heat source, J. Comput. Phys. 228 (2009) 6977e6990.
8
V. Mizonov, N. Yelin / International Journal of Thermal Sciences 97 (2015) 1e8
[11] M.K. Moallemi, R. Visksn, Analysis of melting around a moving heat source, Int. J. Heat Mass Trans. 29 (1986) 1271e1282. [12] H. Berthiaux, V. Mizonov, Applications of Markov chains in particulate process engineering: a review, Can. J. Chem. Eng. 85 (2004) 1143e1168. [13] H. Berthiaux, V. Mizonov, V. Zhukov, Application of the theory of Markov chains to model different processes in particle technology, Powder Technol. 157 (2005) 128e137.
[14] V. Mizonov, H. Berthiaux, P. Arlabosse, D. Djerroud, Application of the theory of Markov chains to model heat and mass transfer between stochastically moving particulate and gas flows, Granul. Matter 10 (2008) 335e340. [15] V. Mizonov, N. Yelin, P. Yakimychev, A cell model to describe and optimize heat and mass transfer in contact heat exchangers, Energy Power Eng. 3 (2011) 144e149.