- Email: [email protected]
Using genetic algorithms for the determination of an heat transfer coefficient in three-phase inverse Stefan problem
Available online at www.sciencedirect.com
International Communications in Heat and Mass Transfer 35 (2008) 149 – 156 www.elsevier.com/locate/ichmt
Using genetic algorithms for the determination of an heat transfer coefficient in three-phase inverse Stefan problem ☆ Damian Słota Institute of Mathematics, Silesian University of Technology, Kaszubska 23, 44-100 Gliwice, Poland Available online 14 September 2007
Abstract In the paper, an example is presented of the application of a genetic algorithm to a design inverse Stefan problem. The problem consists in the reconstruction of the function which describes the heat transfer coefficient, where the positions of phase change moving interfaces are well-known. In numerical calculations, the Tikhonov regularization, a genetic algorithm and a generalized alternating phase truncation method were used. The featured examples of calculations show a very good approximation of the exact solution. © 2007 Elsevier Ltd. All rights reserved. Keywords: Inverse Stefan problem; Tikhonov regularization; Genetic algorithm; Solidification; Heat transfer
1. Introduction The Stefan problem consists in the determination of the temperature distribution within an area and the position of the moving interface of the phase change (freezing front) when the initial condition, boundary conditions and the thermophysical properties of a substance are known [1,2]. The inverse Stefan problem consists in the determination of the initial condition, boundary conditions or the thermophysical properties of a substance [3]. Lack of a portion of the input information is compensated for with extra information about the effects of the input conditions operation. For the inverse Stefan problem, it is assumed that this additional information is the position of the freezing front, its velocity in a normal direction or temperature in selected points of the domain. The problems where the additional information is the position of the moving interface are called design problems. A majority of available studies refer to the one- or two-phase inverse Stefan problem (see for example [4–14]). In papers [15–18], a three-phase design inverse Stefan problem is presented. The problem consists in the reconstruction of the function which describes the heat transfer coefficient, where the positions of the moving interfaces of the phase change are well-known. The authors of paper [15] divide the task into three partial tasks, in turn in a liquid phase, solid phase and in the mushy zone. Next, based on solutions of these tasks, the searched coefficient is determined. The method described in papers [16–18] consists in the minimization of a functional, the value of which is the norm of a
☆
Communicated by W.J. Minkowycz. E-mail address: [email protected].
0735-1933/$ - see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2007.08.010
150
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
Nomenclature A b c d D Dl e J L npop M N pc pm r t T T0 T∞ ⁎ +1 Tk,k Vαp x
coefficients of numerical integration method length of domain (m) specific heat (J·kg− 1 · K− 1) parameter in mutation operator domain of the problem subset of domain D relative percentage error (%) functional latent heat od fusion (J·kg− 1) population size number of control points number of generations crossover probability mutation probability random number time (s) temperature (K) initial temperature (K) ambient temperature (K) phase change temperature (K) set of functions space variable (m)
Greek symbols α heat transfer coefficient (W·m− 2 · K− 1) γ regularization parameter Γl boundary of domain Γk,k + 1 phase change interface ξk,k + 1 phase change interface λ thermal conductivity (W·m− 1 · K− 1) . mass density (kg · m− 3) σ standard deviation σp standard deviation in percentage of the average value τ current generation number
difference between the given positions of the moving interfaces of the phase change and the positions reconstructed from the selected function describing the heat transfer coefficient. In numerical calculations, the Nelder–Mead (downhill simplex) optimization method and the generalized alternating phase truncation method were used. Paper [19] presents a multi-phase design inverse Stefan problem. The solution included in the paper is found in a linear combination form of the functions satisfying the heat conduction equation. The coefficients of the combination are determined by the least square method to minimize the maximal defect in the initial-boundary data. This method facilitates the finding of the temperature distribution, the heat flux or the heat transfer coefficient on the boundary. In this paper, we are going to find a solution of a three-phase design inverse Stefan problem, for which the heat transfer coefficient on one boundary should be so selected that the phase change moving interfaces could take the given positions. The solution will consist in the minimization of the functional whose value is the norm of the difference between the given interfaces' positions and the positions reconstructed for the selected heat transfer coefficient. Stability of the procedure is guaranteed owing to the Tikhonov regularization method [20,21]. For the minimization of
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
151
Fig. 1. Domain of the problem.
the functional, genetic algorithms were used [22–24], whereas the Stefan problem was solved by generalized alternating phase truncation method [18,25,26] (for others methods see for example [27–30]). 2. Formulation of the problem On the boundary of domain D = [0,b] × [0,t⁎] ⊂ ℝ2 three components are distributed (Fig. 1): Γ0 = {(x,0); x ∈ [0,b]}, Γ1 = {(0,t); t ∈ [0,t⁎]}, Γ2 = {(b,t); t ∈ [0,t⁎]}, where initial and boundary conditions are given. Let D1, D2 and D3 be subsets of domain D (D = D1 ∪ D2 ∪ D3) which are occupied by respective phase. Let Γ1,2 will designate the common boundary of domains D1 and D2, whilst Γ2,3 will mean the common boundary of domains D2 and D3. Let us assume that boundary Γk,k + 1 is described by function x = ξk,k + 1(t), k = 1,2. We will look for an approximate solution of the following problem. For the given partial position of interfaces Γ1,2 and Γ2,3, the distribution of temperature Tk in domain Dk (k = 1,2,3) is calculated as well as function α(t) on boundary Γ2, which satisfy the following equations (for k = 1,2,3): ATk 1A ATk kk x c k .k ð x; t Þ ¼ ð x; t Þ ; x Ax At Ax
in Dk ;
ð1Þ
T1 ð x; 0Þ ¼ T0 ;
on C0 ;
ð2Þ
ATk ð x; t Þ ¼ 0; Ax
on C1 ;
ð3Þ
on C2 ;
ð4Þ
on Ck;kþ1 ;
ð5Þ
−kk
ATk ð x; t Þ ¼ aðt ÞðTk ð x; t Þ Tl Þ; Ax
T ; Tk nk;kþ1 ðt Þ; t ¼ Tkþ1 nk;kþ1 ðt Þ; t ¼ Tk;kþ1
Lk;kþ1 .kþ1
dnk;kþ1 ðt Þ ATk ð x; t Þ ¼ kk Ax dt
j
Ck;kþ1
þ kkþ1
ATkþ1 ð x; t Þ Ax
j
Ck;kþ1
;
ð6Þ
152
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
where ck are the specific heat in respective phase, λk are the thermal conductivity, ϱk are the mass density, α is the heat ⁎ + 1 are the phase change transfer coefficient, T0 is the initial temperature, T∞ is the ambient temperature, Tk,k temperature, Lk,k + 1 are the latent heat of fusion, and x and t refer to spatial location and time, respectively. We will look for the α(t) function in the form: 8 a1 > > > > > < a2 aðt Þ ¼ v > > > an1 > > : an
for t V t1 for t a ðt1 ; t2 ð7Þ for t a ðtn2 ;tn1 for t N tn1 ;
where 0 b t1 b t2 b … b tn − 1 b t⁎. Let Vαp mean a set of all functions in the form (Eq. (7)), where αi ∈ [αil, αiu], αil b αiu, αil, αiu ∈ ℝ, for i = 1, 2,…, n. For the given function α(t) ∈ Vαp, problems (1)–(6) becomes a direct three-phase Stefan problem, whose solution enables finding the positions of the moving interfaces corresponding to the α(t) function. We will assume in the calculations, that the positions of the moving interfaces are given in selected points only (hereafter called the “control points”). Using the found interfaces' positions, ξk,k + 1(t), and the given positions ξ*k,k + 1(t), k = 1, 2, we can build a functional which will specify the error of an approximated solution: J ðaðt ÞÞ ¼
2 X k¼1
jj n
k;kþ1 ðt Þ
n4k;kþ1 ðt Þ
jj
2
jj
þ g að t Þ
jj
2
ð8Þ
where γ is the regularization parameter, determined based on the discrepancy principle [20,21]. The above norms designate the norms in a space of square-integrable functions in the interval (0, t⁎) (L2(0, t⁎)). When the exact position of the interface is given only in selected points, the first norm present in functional (8) will be calculated from the following dependence:
jj
nk;kþ1 ðt Þ n4k;kþ1 ðt Þ
jj
2
¼
M 2 X Ai nk;kþ1;i n4k;kþ1;i
ð9Þ
i¼1
where Ai are coefficients dependent on the chosen numerical integration method, M is the number of control points, and ⁎ ξ⁎ k,k + 1;i = ξk,k + 1 (ti) and ξk,k + 1;i = ξk,k + 1(ti) are the given and calculated points respectively, describing the interfaces positions. 3. Genetic algorithm For the minimization of the functional genetic algorithms were used. In which for the representation of the vector of decision variables (α1,…,αn), a chromosome was used in the form of a vector of the real numbers (real number representation) [22–24]. A tournament selection was applied in the algorithm. The selection is carried out so that two chromosomes are drawn and the one with better fitness goes to a new generation. There are as many draws as individuals that the new generation is supposed to include. An elitist model was also applied in the algorithm, where the best individual of the previous generation is saved and, if all individuals in the current generation are worse, the worst of them is replaced with the saved best individual from the previous population. As the crossover operator, arithmetical crossover was applied, where as a result of crossing of chromosomes α1 = (α11,…, α1n) and α2 = (α12,…, αn2), their linear combinations are obtained: α1 V ¼ rα1 þ ð1 rÞα2 α2 V ¼ rα2 þ ð1 rÞα1 where parameter r is a random number with a uniform distribution from the domain [0,1].
ð10Þ ð11Þ
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
153
In the calculations, a nonuniform mutation operator was used as well. During mutation, the αi gene is transformed according to the equation:
ai þ D s; aui ai V ai ¼ ð12Þ ai D s; ai ali and a decision is taken at random about which from the above formulas should be applied. Function Δ(τ, x) was assumed in the form: s Dðs; xÞ ¼ x 1 rð1N Þd ð13Þ where r is a random number with a uniform distribution from the domain [0, 1], τ is the current generation number, N is the maximum number of generations and d is a constant parameter (in the calculations, d = 2 was assumed). 4. Calculations It was assumed in the calculations that: b = 0.06 [m], λ1 = 54 [W·m− 1·K− 1], λ2 = 30 [W·m− 1·K− 1], λ3 = 42 [W·m− 1·K− 1], c1 = 840 [J·kg− 1·K− 1], c2 = 670 [J·kg− 1·K− 1], c3 = 755 [J·kg− 1·K− 1], .1 = 7000 [kg·m− 3], .2 = 7500 ⁎ = 1773 [K] and T2,3 ⁎ = 1718 [K], The [kg·m− 3], .3 = 7250 [kg·m− 3], L1,2 = 217600 [J·kg− 1], L2,3 = 54400 [J·kg− 1], T1,2 ambient temperature equal T∞ = 303 [K] and initial temperature is equal T0 = 1803 [K]. The exact value of the heat transfer coefficient amounts to: 8 1460 > > > > > < 872 að zÞ ¼ 695 > > > > 515 > : 250
for t V t1 for t a ðt1 ; t2 for t a ðt2 ; t3
ð14Þ
for t a ðt3 ; t4 for t N t4
where t1 = 10 [s], t2 = 35 [s], t3 = 55 [s], t4 = 85 [s]. The parameters used for the genetic algorithm are as follows: population size npop = 70, number of generations N = 1000, crossover probability pc = 0.7 and mutation probability pm = 0.1. In the generalized alternating phase truncation method, the finite-difference method was used, the calculations having been made on a grid of discretization intervals equal Δt = 1/10 and Δx = b/500. A reasonable change of the grid density did not significantly affect the results obtained. Calculations were carried out for five different initial settings of a pseudorandom numbers' generator. In the example, the influence was studied of the number of control points, i.e. the number of points in which the moving interfaces positions are known, on accuracy of the heat transfer coefficient determination. The calculations were made with control of the moving interfaces positions, conducted every second, every two and every five seconds. This corresponds to a situation where the sum occurring in the functional (8) consists of 243, 121 or 49 components, respectively. Also, the influence was investigated of the measuring errors with which the predefined positions of the moving interfaces are burdened on accuracy of the calculations. In order to examine the influence of measuring errors on the accuracy of the results obtained, thus to examine the stability of the algorithm, the solution of the direct Stefan problem was disturbed with a pseudorandom error of normal distribution. Next, the so obtained positions of the moving interfaces were used as input information for the inverse Stefan problem. Tables 1–3 show the results of the heat transfer coefficient calculations and average values of the relative percentage errors with which parameters αi were reconstructed. These tables juxtapose the values of standard deviation obtained for calculations made at different initial settings of the pseudorandom numbers generator. The last column of the tables includes the standard deviation values in percentage of the average value. As shown in the results, the presented algorithm is stable in terms of the input data errors. Each time when the input data were burdened with errors, the error of the heat transfer coefficient reconstruction did not exceed the initial error, whereas in the case of undisturbed input data, the errors of results are insignificant. In this case, the scatter or results is the smallest and it grows as the input data errors increase. Also, a (reasonable) decrease of the number of control points
154
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
Table 1 Results of the boundary condition reconstruction for interfaces control every one second Per. (%)
α
Error (%)
σ
σp (%)
0
α1 = 1460.153 α2 = 871.841 α3 = 694.432 α4 = 515.793 α5 = 249.856 α1 = 1458.910 α2 = 871.824 α3 = 693.092 α4 = 518.227 α5 = 249.360 α1 = 1453.742 α2 = 875.902 α3 = 691.852 α4 = 516.451 α5 = 249.603
e1 = 0.010479 e2 = 0.018234 e3 = 0.081727 e4 = 0.153981 e5 = 0.057600 e1 = 0.074658 e2 = 0.020183 e3 = 0.274532 e4 = 0.626602 e5 = 0.256000 e1 = 0.428630 e2 = 0.447477 e3 = 0.452950 e4 = 0.281748 e5 = 0.158800
σ1 = 0.307173 σ2 = 1.407981 σ3 = 3.679839 σ4 = 1.691143 σ5 = 0.165965 σ1 = 0.868378 σ2 = 2.831365 σ3 = 8.139637 σ4 = 4.352580 σ5 = 0.492819 σ1 = 5.524236 σ2 = 3.956776 σ3 = 6.999785 σ4 = 4.658699 σ5 = 0.803302
σp1 = 0.021037 σp2 = 0.161495 σp3 = 0.529906 σp4 = 0.327872 σp5 = 0.066424 σp1 = 0.059522 σp2 = 0.324763 σp3 = 1.174396 σp4 = 0.839899 σp5 = 0.197633 σp1 = 0.380001 σp2 = 0.451737 σp3 = 1.011746 σp4 = 0.902061 σp5 = 0.321831
1
2
Table 2 Results of the boundary condition reconstruction for interfaces control every two seconds Per. (%)
α
Error (%)
σ
σp (%)
0
α1 = 1459.111 α2 = 871.505 α3 = 695.268 α4 = 515.492 α5 = 250.057 α1 = 1458.074 α2 = 873.198 α3 = 691.722 α4 = 517.023 α5 = 249.869 α1 = 1446.634 α2 = 877.140 α3 = 685.273 α4 = 523.375 α5 = 248.937
e1 = 0.060890 e2 = 0.056766 e3 = 0.038561 e4 = 0.095534 e5 = 0.022800 e1 = 0.131918 e2 = 0.137385 e3 = 0.471655 e4 = 0.392816 e5 = 0.052400 e1 = 0.915479 e2 = 0.589450 e3 = 1.399568 e4 = 1.626214 e5 = 0.425200
σ1 = 1.115842 σ2 = 1.919349 σ3 = 4.813095 σ4 = 2.147596 σ5 = 0.134666 σ1 = 0.723096 σ2 = 1.507263 σ3 = 3.183905 σ4 = 1.092505 σ5 = 0.114190 σ1 = 15.533968 σ2 = 10.890536 σ3 = 11.505990 σ4 = 5.3203841 σ5 = 0.8257204
σp1 = 0.076474 σp2 = 0.220234 σp3 = 0.692265 σp4 = 0.416611 σp5 = 0.053854 σp1 = 0.049593 σp2 = 0.172614 σp3 = 0.460286 σp4 = 0.211307 σp5 = 0.045700 σp1 = 1.073801 σp2 = 1.241596 σp3 = 1.679038 σp4 = 1.016553 σp5 = 0.331699
1
2
Table 3 Results of the boundary condition reconstruction for interfaces control every five seconds Per. (%)
α
Error (%)
σ
σp (%)
0
α1 = 1461.131 α2 = 871.707 α3 = 695.265 α4 = 514.662 α5 = 250.231 α1 = 1461.308 α2 = 870.468 α3 = 697.673 α4 = 513.494 α5 = 249.841 α1 = 1448.023 α2 = 876.608 α3 = 694.550 α4 = 513.718 α5 = 250.562
e1 = 0.077466 e2 = 0.033601 e3 = 0.038129 e4 = 0.065631 e5 = 0.092400 e1 = 0.089589 e2 = 0.175688 e3 = 0.384604 e4 = 0.292427 e5 = 0.063600 e1 = 0.820342 e2 = 0.528440 e3 = 0.064748 e4 = 0.248932 e5 = 0.224800
σ1 = 1.446537 σ2 = 1.203757 σ3 = 2.110147 σ4 = 0.898082 σ5 = 0.221733 σ1 = 4.602703 σ2 = 4.822104 σ3 = 3.404690 σ4 = 1.559828 σ5 = 0.763081 σ1 = 16.269957 σ2 = 11.109687 σ3 = 14.676119 σ4 = 7.4259390 σ5 = 1.4771168
σp1 = 0.099001 σp2 = 0.138092 σp3 = 0.303503 σp4 = 0.174500 σp5 = 0.088611 σp1 = 0.314971 σp2 = 0.553967 σp3 = 0.488007 σp4 = 0.303768 σp5 = 0.305427 σp1 = 1.123598 σp2 = 1.267350 σp3 = 2.113041 σp4 = 1.445528 σp5 = 0.589521
1
2
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
155
Fig. 2. Exact (solid lines) and reconstructed (dot lines) positions of the moving interfaces for perturbation equal to 2% and interfaces control every five seconds.
does not result in a loss of the algorithm's stability. Fig. 2 presents the accurate and reconstructed positions of the moving interfaces, obtained for input data disturbed with a 2% error for interfaces control carried out every five seconds. In the remaining cases, the interfaces were reconstructed equally well. Fig. 3 shows a comparison of the heat transfer coefficient reconstruction errors obtained via a procedure which uses genetic algorithms and via a procedure using the Nelder–Mead method [16] (for the Nelder–Mead method see for example [31]). As can be seen, the application of genetic algorithms has considerably improved the results obtained. The same effect (although the differences were much smaller) was noted in the case of a two-phase inverse Stefan problem [10]. In the case of the Nelder–Mead method, every subsequent actuation of the method with different initial points yielded a larger scatter of results than the subsequent actuations of the genetic algorithm (see [10]). 5. Conclusion This paper discussed the identification of the heat transfer coefficient in the inverse three-phase Stefan problem. The problem consists in the reconstruction of the function which describes the heat transfer coefficient, where the position of the moving interfaces of the phase change are well-known (inverse design Stefan problem). In numerical calculations, a genetic algorithm, the Tikhonov regularization and a generalized alternating phase truncation method were used. The calculations made show stability of the algorithm in terms of the number of control points and the input data errors. The application of genetic algorithms enhanced the accuracy of the results obtained, in relation to the results obtained when using the Nelder–Mead method. What is also important, is the small scatter of the results obtained during calculations for different initial settings of the pseudorandom numbers generator.
Fig. 3. Comparison of average (a) and maximum (b) errors of the heat transfer coefficient reconstruction, obtained with the use of genetic algorithms (GA) and the Nelder–Mead method (NM) for a different number of control points.
156
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
An example has been presented in the paper of the heat transfer coefficient selection in a few specified zones. The task can be transformed without any problem so that the amount and length of zones would be selected at a defined (and variable) value of the heat transfer coefficient. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]
L.I. Rubinstein, The Stefan Problem, AMS, Providence, 1971. A.M. Meirmanov, The Stefan Problem, Walter de Gruyter, Berlin, 1992. N.L. Goldman, Inverse Stefan Problem, Kluwer, Dordrecht, 1997. N. Al-Khalidy, Application of optimization methods for solving inverse phase-change problems, Numer. Heat Transf., B Fundam. 31 (1997) 477–497. D.D. Ang, A.P.N. Dinh, D.N. Thanh, Regularization of an inverse two-phase Stefan problem, Nonlinear Anal. 34 (1998) 719–731. R. Grzymkowski, D. Słota, One-phase inverse Stefan problems solved by Adomian decomposition method, Comput. Math. Appl. 51 (2006) 33–40. P. Jochum, The numerical solution of the inverse Stefan problem, Numer. Math. 34 (1980) 411–429. J. Liu, B. Guerrier, A comparative study of domain embedding methods for regularized solutions of inverse Stefan problems, Int. J. Numer. Methods Eng. 40 (1997) 3579–3600. D. Słota, Influence of the mutation operator on the solution of an inverse Stefan problem by genetic algorithms, in: V.N. Alexandrov, et al., (Eds.), Computational Science — ICCS 2006, LNCS 3991, Springer-Verlag, Berlin, 2006, pp. 786–789. D. Słota, Solving the inverse Stefan design problem using genetic algorithms, Inverse Probl. Sci. Eng. (in review). V.R. Voller, Enthalpy method for inverse Stefan problems, Numer. Heat Transf., B Fundam. 21 (1992) 41–55. N. Zabaras, Y. Ruan, O. Richmond, Design of two-dimensional Stefan processes with desired freezing front motions, Numer. Heat Transf., B Fundam. 21 (1992) 307–325. N. Zabaras, S. Kang, On the solution of an ill-posed design solidification problem using minimization techniques in finite- and infinitedimensional function space, Int. J. Numer. Methods Eng. 36 (1993) 3973–3990. N. Zabaras, K. Yuan, Dynamic programming approach to the inverse Stefan design problem, Numer. Heat Transf., B Fundam. 26 (1994) 97–104. M. Slodička, H.D. Schepper, Determination of the heat-transfer coefficient during soldification of alloys, Comput. Methods Appl. Mech. Eng. 194 (2005) 491–498. R. Grzymkowski, D. Słota, Determination of unknown thermal coefficient in continuous casting, in: H.G. Fritz (Ed.), Proc. of the 3rd ESAFORM Conference on Material Forming, Universität Stuttgart, Stuttgart, 2000, pp. X11–X14. R. Grzymkowski, D. Słota, Numerical calculations of the heat-transfer coefficient during solidification of alloys, in: B. Sarler, C.A. Brebbia (Eds.), Moving Boundaries VI, Wit Press, Southampton, 2001, pp. 41–50. R. Grzymkowski, D. Słota, Numerical method for multi-phase inverse Stefan design problems, Arch. Metall. Mater. 51 (2006) 161–172. R. Grzymkowski, D. Słota, Multi-phase inverse Stefan problems solved by approximation method, in: R. Wyrzykowski, et al., (Eds.), Parallel Processing and Applied Mathematic, LNCS 2328, Springer-Verlag, Berlin, 2002, pp. 679–686. K. Kurpisz, A.J. Nowak, Inverse Thermal Problems, Comput. Mech. Publ., Southampton, 1995. A.N. Tikhonov, V.Y. Arsenin, Solution of Ill-Posed Problems, Wiley & Sons, New York, 1977. D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Reading, 1989. Z. Michalewicz, Genetic Algorithms + Data Structures = Evolution Programs, Springer-Verlag, Berlin, 1996. A. Osyczka, Evolutionary Algorithms for Single and Multicriteria Design Optimization, Physica-Verlag, Heidelberg, 2002. A. Kapusta, B. Mochnacki, The analysis of heat transfer processes in the cylindrical radial continuous casting volume, Bull. Pol. Acad. Sci., Tech. Sci. 36 (1988) 309–320. E. Majchrzak, B. Mochnacki, Application of the BEM in the thermal theory of foundry, Eng. Anal. Bound. Elem. 16 (1995) 99–121. J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1996. B. Furenes, B. Lie, Using event location in finite-difference methods for phase-change problems, Numer. Heat Transf., B Fundam. 50 (2006) 143–155. Rizwan-uddin, One-dimensional phase change with periodic boundary conditions, Numer. Heat Transf., A Appl. 35 (1999) 361–372. M. Zerroukat, C. R. Chatwin, Computational Moving Boundary Problems, Research Studies Press, Taunton, 1994. J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313.
International Communications in Heat and Mass Transfer 35 (2008) 149 – 156 www.elsevier.com/locate/ichmt
Using genetic algorithms for the determination of an heat transfer coefficient in three-phase inverse Stefan problem ☆ Damian Słota Institute of Mathematics, Silesian University of Technology, Kaszubska 23, 44-100 Gliwice, Poland Available online 14 September 2007
Abstract In the paper, an example is presented of the application of a genetic algorithm to a design inverse Stefan problem. The problem consists in the reconstruction of the function which describes the heat transfer coefficient, where the positions of phase change moving interfaces are well-known. In numerical calculations, the Tikhonov regularization, a genetic algorithm and a generalized alternating phase truncation method were used. The featured examples of calculations show a very good approximation of the exact solution. © 2007 Elsevier Ltd. All rights reserved. Keywords: Inverse Stefan problem; Tikhonov regularization; Genetic algorithm; Solidification; Heat transfer
1. Introduction The Stefan problem consists in the determination of the temperature distribution within an area and the position of the moving interface of the phase change (freezing front) when the initial condition, boundary conditions and the thermophysical properties of a substance are known [1,2]. The inverse Stefan problem consists in the determination of the initial condition, boundary conditions or the thermophysical properties of a substance [3]. Lack of a portion of the input information is compensated for with extra information about the effects of the input conditions operation. For the inverse Stefan problem, it is assumed that this additional information is the position of the freezing front, its velocity in a normal direction or temperature in selected points of the domain. The problems where the additional information is the position of the moving interface are called design problems. A majority of available studies refer to the one- or two-phase inverse Stefan problem (see for example [4–14]). In papers [15–18], a three-phase design inverse Stefan problem is presented. The problem consists in the reconstruction of the function which describes the heat transfer coefficient, where the positions of the moving interfaces of the phase change are well-known. The authors of paper [15] divide the task into three partial tasks, in turn in a liquid phase, solid phase and in the mushy zone. Next, based on solutions of these tasks, the searched coefficient is determined. The method described in papers [16–18] consists in the minimization of a functional, the value of which is the norm of a
☆
Communicated by W.J. Minkowycz. E-mail address: [email protected].
0735-1933/$ - see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2007.08.010
150
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
Nomenclature A b c d D Dl e J L npop M N pc pm r t T T0 T∞ ⁎ +1 Tk,k Vαp x
coefficients of numerical integration method length of domain (m) specific heat (J·kg− 1 · K− 1) parameter in mutation operator domain of the problem subset of domain D relative percentage error (%) functional latent heat od fusion (J·kg− 1) population size number of control points number of generations crossover probability mutation probability random number time (s) temperature (K) initial temperature (K) ambient temperature (K) phase change temperature (K) set of functions space variable (m)
Greek symbols α heat transfer coefficient (W·m− 2 · K− 1) γ regularization parameter Γl boundary of domain Γk,k + 1 phase change interface ξk,k + 1 phase change interface λ thermal conductivity (W·m− 1 · K− 1) . mass density (kg · m− 3) σ standard deviation σp standard deviation in percentage of the average value τ current generation number
difference between the given positions of the moving interfaces of the phase change and the positions reconstructed from the selected function describing the heat transfer coefficient. In numerical calculations, the Nelder–Mead (downhill simplex) optimization method and the generalized alternating phase truncation method were used. Paper [19] presents a multi-phase design inverse Stefan problem. The solution included in the paper is found in a linear combination form of the functions satisfying the heat conduction equation. The coefficients of the combination are determined by the least square method to minimize the maximal defect in the initial-boundary data. This method facilitates the finding of the temperature distribution, the heat flux or the heat transfer coefficient on the boundary. In this paper, we are going to find a solution of a three-phase design inverse Stefan problem, for which the heat transfer coefficient on one boundary should be so selected that the phase change moving interfaces could take the given positions. The solution will consist in the minimization of the functional whose value is the norm of the difference between the given interfaces' positions and the positions reconstructed for the selected heat transfer coefficient. Stability of the procedure is guaranteed owing to the Tikhonov regularization method [20,21]. For the minimization of
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
151
Fig. 1. Domain of the problem.
the functional, genetic algorithms were used [22–24], whereas the Stefan problem was solved by generalized alternating phase truncation method [18,25,26] (for others methods see for example [27–30]). 2. Formulation of the problem On the boundary of domain D = [0,b] × [0,t⁎] ⊂ ℝ2 three components are distributed (Fig. 1): Γ0 = {(x,0); x ∈ [0,b]}, Γ1 = {(0,t); t ∈ [0,t⁎]}, Γ2 = {(b,t); t ∈ [0,t⁎]}, where initial and boundary conditions are given. Let D1, D2 and D3 be subsets of domain D (D = D1 ∪ D2 ∪ D3) which are occupied by respective phase. Let Γ1,2 will designate the common boundary of domains D1 and D2, whilst Γ2,3 will mean the common boundary of domains D2 and D3. Let us assume that boundary Γk,k + 1 is described by function x = ξk,k + 1(t), k = 1,2. We will look for an approximate solution of the following problem. For the given partial position of interfaces Γ1,2 and Γ2,3, the distribution of temperature Tk in domain Dk (k = 1,2,3) is calculated as well as function α(t) on boundary Γ2, which satisfy the following equations (for k = 1,2,3): ATk 1A ATk kk x c k .k ð x; t Þ ¼ ð x; t Þ ; x Ax At Ax
in Dk ;
ð1Þ
T1 ð x; 0Þ ¼ T0 ;
on C0 ;
ð2Þ
ATk ð x; t Þ ¼ 0; Ax
on C1 ;
ð3Þ
on C2 ;
ð4Þ
on Ck;kþ1 ;
ð5Þ
−kk
ATk ð x; t Þ ¼ aðt ÞðTk ð x; t Þ Tl Þ; Ax
T ; Tk nk;kþ1 ðt Þ; t ¼ Tkþ1 nk;kþ1 ðt Þ; t ¼ Tk;kþ1
Lk;kþ1 .kþ1
dnk;kþ1 ðt Þ ATk ð x; t Þ ¼ kk Ax dt
j
Ck;kþ1
þ kkþ1
ATkþ1 ð x; t Þ Ax
j
Ck;kþ1
;
ð6Þ
152
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
where ck are the specific heat in respective phase, λk are the thermal conductivity, ϱk are the mass density, α is the heat ⁎ + 1 are the phase change transfer coefficient, T0 is the initial temperature, T∞ is the ambient temperature, Tk,k temperature, Lk,k + 1 are the latent heat of fusion, and x and t refer to spatial location and time, respectively. We will look for the α(t) function in the form: 8 a1 > > > > > < a2 aðt Þ ¼ v > > > an1 > > : an
for t V t1 for t a ðt1 ; t2 ð7Þ for t a ðtn2 ;tn1 for t N tn1 ;
where 0 b t1 b t2 b … b tn − 1 b t⁎. Let Vαp mean a set of all functions in the form (Eq. (7)), where αi ∈ [αil, αiu], αil b αiu, αil, αiu ∈ ℝ, for i = 1, 2,…, n. For the given function α(t) ∈ Vαp, problems (1)–(6) becomes a direct three-phase Stefan problem, whose solution enables finding the positions of the moving interfaces corresponding to the α(t) function. We will assume in the calculations, that the positions of the moving interfaces are given in selected points only (hereafter called the “control points”). Using the found interfaces' positions, ξk,k + 1(t), and the given positions ξ*k,k + 1(t), k = 1, 2, we can build a functional which will specify the error of an approximated solution: J ðaðt ÞÞ ¼
2 X k¼1
jj n
k;kþ1 ðt Þ
n4k;kþ1 ðt Þ
jj
2
jj
þ g að t Þ
jj
2
ð8Þ
where γ is the regularization parameter, determined based on the discrepancy principle [20,21]. The above norms designate the norms in a space of square-integrable functions in the interval (0, t⁎) (L2(0, t⁎)). When the exact position of the interface is given only in selected points, the first norm present in functional (8) will be calculated from the following dependence:
jj
nk;kþ1 ðt Þ n4k;kþ1 ðt Þ
jj
2
¼
M 2 X Ai nk;kþ1;i n4k;kþ1;i
ð9Þ
i¼1
where Ai are coefficients dependent on the chosen numerical integration method, M is the number of control points, and ⁎ ξ⁎ k,k + 1;i = ξk,k + 1 (ti) and ξk,k + 1;i = ξk,k + 1(ti) are the given and calculated points respectively, describing the interfaces positions. 3. Genetic algorithm For the minimization of the functional genetic algorithms were used. In which for the representation of the vector of decision variables (α1,…,αn), a chromosome was used in the form of a vector of the real numbers (real number representation) [22–24]. A tournament selection was applied in the algorithm. The selection is carried out so that two chromosomes are drawn and the one with better fitness goes to a new generation. There are as many draws as individuals that the new generation is supposed to include. An elitist model was also applied in the algorithm, where the best individual of the previous generation is saved and, if all individuals in the current generation are worse, the worst of them is replaced with the saved best individual from the previous population. As the crossover operator, arithmetical crossover was applied, where as a result of crossing of chromosomes α1 = (α11,…, α1n) and α2 = (α12,…, αn2), their linear combinations are obtained: α1 V ¼ rα1 þ ð1 rÞα2 α2 V ¼ rα2 þ ð1 rÞα1 where parameter r is a random number with a uniform distribution from the domain [0,1].
ð10Þ ð11Þ
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
153
In the calculations, a nonuniform mutation operator was used as well. During mutation, the αi gene is transformed according to the equation:
ai þ D s; aui ai V ai ¼ ð12Þ ai D s; ai ali and a decision is taken at random about which from the above formulas should be applied. Function Δ(τ, x) was assumed in the form: s Dðs; xÞ ¼ x 1 rð1N Þd ð13Þ where r is a random number with a uniform distribution from the domain [0, 1], τ is the current generation number, N is the maximum number of generations and d is a constant parameter (in the calculations, d = 2 was assumed). 4. Calculations It was assumed in the calculations that: b = 0.06 [m], λ1 = 54 [W·m− 1·K− 1], λ2 = 30 [W·m− 1·K− 1], λ3 = 42 [W·m− 1·K− 1], c1 = 840 [J·kg− 1·K− 1], c2 = 670 [J·kg− 1·K− 1], c3 = 755 [J·kg− 1·K− 1], .1 = 7000 [kg·m− 3], .2 = 7500 ⁎ = 1773 [K] and T2,3 ⁎ = 1718 [K], The [kg·m− 3], .3 = 7250 [kg·m− 3], L1,2 = 217600 [J·kg− 1], L2,3 = 54400 [J·kg− 1], T1,2 ambient temperature equal T∞ = 303 [K] and initial temperature is equal T0 = 1803 [K]. The exact value of the heat transfer coefficient amounts to: 8 1460 > > > > > < 872 að zÞ ¼ 695 > > > > 515 > : 250
for t V t1 for t a ðt1 ; t2 for t a ðt2 ; t3
ð14Þ
for t a ðt3 ; t4 for t N t4
where t1 = 10 [s], t2 = 35 [s], t3 = 55 [s], t4 = 85 [s]. The parameters used for the genetic algorithm are as follows: population size npop = 70, number of generations N = 1000, crossover probability pc = 0.7 and mutation probability pm = 0.1. In the generalized alternating phase truncation method, the finite-difference method was used, the calculations having been made on a grid of discretization intervals equal Δt = 1/10 and Δx = b/500. A reasonable change of the grid density did not significantly affect the results obtained. Calculations were carried out for five different initial settings of a pseudorandom numbers' generator. In the example, the influence was studied of the number of control points, i.e. the number of points in which the moving interfaces positions are known, on accuracy of the heat transfer coefficient determination. The calculations were made with control of the moving interfaces positions, conducted every second, every two and every five seconds. This corresponds to a situation where the sum occurring in the functional (8) consists of 243, 121 or 49 components, respectively. Also, the influence was investigated of the measuring errors with which the predefined positions of the moving interfaces are burdened on accuracy of the calculations. In order to examine the influence of measuring errors on the accuracy of the results obtained, thus to examine the stability of the algorithm, the solution of the direct Stefan problem was disturbed with a pseudorandom error of normal distribution. Next, the so obtained positions of the moving interfaces were used as input information for the inverse Stefan problem. Tables 1–3 show the results of the heat transfer coefficient calculations and average values of the relative percentage errors with which parameters αi were reconstructed. These tables juxtapose the values of standard deviation obtained for calculations made at different initial settings of the pseudorandom numbers generator. The last column of the tables includes the standard deviation values in percentage of the average value. As shown in the results, the presented algorithm is stable in terms of the input data errors. Each time when the input data were burdened with errors, the error of the heat transfer coefficient reconstruction did not exceed the initial error, whereas in the case of undisturbed input data, the errors of results are insignificant. In this case, the scatter or results is the smallest and it grows as the input data errors increase. Also, a (reasonable) decrease of the number of control points
154
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
Table 1 Results of the boundary condition reconstruction for interfaces control every one second Per. (%)
α
Error (%)
σ
σp (%)
0
α1 = 1460.153 α2 = 871.841 α3 = 694.432 α4 = 515.793 α5 = 249.856 α1 = 1458.910 α2 = 871.824 α3 = 693.092 α4 = 518.227 α5 = 249.360 α1 = 1453.742 α2 = 875.902 α3 = 691.852 α4 = 516.451 α5 = 249.603
e1 = 0.010479 e2 = 0.018234 e3 = 0.081727 e4 = 0.153981 e5 = 0.057600 e1 = 0.074658 e2 = 0.020183 e3 = 0.274532 e4 = 0.626602 e5 = 0.256000 e1 = 0.428630 e2 = 0.447477 e3 = 0.452950 e4 = 0.281748 e5 = 0.158800
σ1 = 0.307173 σ2 = 1.407981 σ3 = 3.679839 σ4 = 1.691143 σ5 = 0.165965 σ1 = 0.868378 σ2 = 2.831365 σ3 = 8.139637 σ4 = 4.352580 σ5 = 0.492819 σ1 = 5.524236 σ2 = 3.956776 σ3 = 6.999785 σ4 = 4.658699 σ5 = 0.803302
σp1 = 0.021037 σp2 = 0.161495 σp3 = 0.529906 σp4 = 0.327872 σp5 = 0.066424 σp1 = 0.059522 σp2 = 0.324763 σp3 = 1.174396 σp4 = 0.839899 σp5 = 0.197633 σp1 = 0.380001 σp2 = 0.451737 σp3 = 1.011746 σp4 = 0.902061 σp5 = 0.321831
1
2
Table 2 Results of the boundary condition reconstruction for interfaces control every two seconds Per. (%)
α
Error (%)
σ
σp (%)
0
α1 = 1459.111 α2 = 871.505 α3 = 695.268 α4 = 515.492 α5 = 250.057 α1 = 1458.074 α2 = 873.198 α3 = 691.722 α4 = 517.023 α5 = 249.869 α1 = 1446.634 α2 = 877.140 α3 = 685.273 α4 = 523.375 α5 = 248.937
e1 = 0.060890 e2 = 0.056766 e3 = 0.038561 e4 = 0.095534 e5 = 0.022800 e1 = 0.131918 e2 = 0.137385 e3 = 0.471655 e4 = 0.392816 e5 = 0.052400 e1 = 0.915479 e2 = 0.589450 e3 = 1.399568 e4 = 1.626214 e5 = 0.425200
σ1 = 1.115842 σ2 = 1.919349 σ3 = 4.813095 σ4 = 2.147596 σ5 = 0.134666 σ1 = 0.723096 σ2 = 1.507263 σ3 = 3.183905 σ4 = 1.092505 σ5 = 0.114190 σ1 = 15.533968 σ2 = 10.890536 σ3 = 11.505990 σ4 = 5.3203841 σ5 = 0.8257204
σp1 = 0.076474 σp2 = 0.220234 σp3 = 0.692265 σp4 = 0.416611 σp5 = 0.053854 σp1 = 0.049593 σp2 = 0.172614 σp3 = 0.460286 σp4 = 0.211307 σp5 = 0.045700 σp1 = 1.073801 σp2 = 1.241596 σp3 = 1.679038 σp4 = 1.016553 σp5 = 0.331699
1
2
Table 3 Results of the boundary condition reconstruction for interfaces control every five seconds Per. (%)
α
Error (%)
σ
σp (%)
0
α1 = 1461.131 α2 = 871.707 α3 = 695.265 α4 = 514.662 α5 = 250.231 α1 = 1461.308 α2 = 870.468 α3 = 697.673 α4 = 513.494 α5 = 249.841 α1 = 1448.023 α2 = 876.608 α3 = 694.550 α4 = 513.718 α5 = 250.562
e1 = 0.077466 e2 = 0.033601 e3 = 0.038129 e4 = 0.065631 e5 = 0.092400 e1 = 0.089589 e2 = 0.175688 e3 = 0.384604 e4 = 0.292427 e5 = 0.063600 e1 = 0.820342 e2 = 0.528440 e3 = 0.064748 e4 = 0.248932 e5 = 0.224800
σ1 = 1.446537 σ2 = 1.203757 σ3 = 2.110147 σ4 = 0.898082 σ5 = 0.221733 σ1 = 4.602703 σ2 = 4.822104 σ3 = 3.404690 σ4 = 1.559828 σ5 = 0.763081 σ1 = 16.269957 σ2 = 11.109687 σ3 = 14.676119 σ4 = 7.4259390 σ5 = 1.4771168
σp1 = 0.099001 σp2 = 0.138092 σp3 = 0.303503 σp4 = 0.174500 σp5 = 0.088611 σp1 = 0.314971 σp2 = 0.553967 σp3 = 0.488007 σp4 = 0.303768 σp5 = 0.305427 σp1 = 1.123598 σp2 = 1.267350 σp3 = 2.113041 σp4 = 1.445528 σp5 = 0.589521
1
2
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
155
Fig. 2. Exact (solid lines) and reconstructed (dot lines) positions of the moving interfaces for perturbation equal to 2% and interfaces control every five seconds.
does not result in a loss of the algorithm's stability. Fig. 2 presents the accurate and reconstructed positions of the moving interfaces, obtained for input data disturbed with a 2% error for interfaces control carried out every five seconds. In the remaining cases, the interfaces were reconstructed equally well. Fig. 3 shows a comparison of the heat transfer coefficient reconstruction errors obtained via a procedure which uses genetic algorithms and via a procedure using the Nelder–Mead method [16] (for the Nelder–Mead method see for example [31]). As can be seen, the application of genetic algorithms has considerably improved the results obtained. The same effect (although the differences were much smaller) was noted in the case of a two-phase inverse Stefan problem [10]. In the case of the Nelder–Mead method, every subsequent actuation of the method with different initial points yielded a larger scatter of results than the subsequent actuations of the genetic algorithm (see [10]). 5. Conclusion This paper discussed the identification of the heat transfer coefficient in the inverse three-phase Stefan problem. The problem consists in the reconstruction of the function which describes the heat transfer coefficient, where the position of the moving interfaces of the phase change are well-known (inverse design Stefan problem). In numerical calculations, a genetic algorithm, the Tikhonov regularization and a generalized alternating phase truncation method were used. The calculations made show stability of the algorithm in terms of the number of control points and the input data errors. The application of genetic algorithms enhanced the accuracy of the results obtained, in relation to the results obtained when using the Nelder–Mead method. What is also important, is the small scatter of the results obtained during calculations for different initial settings of the pseudorandom numbers generator.
Fig. 3. Comparison of average (a) and maximum (b) errors of the heat transfer coefficient reconstruction, obtained with the use of genetic algorithms (GA) and the Nelder–Mead method (NM) for a different number of control points.
156
D. Słota / International Communications in Heat and Mass Transfer 35 (2008) 149–156
An example has been presented in the paper of the heat transfer coefficient selection in a few specified zones. The task can be transformed without any problem so that the amount and length of zones would be selected at a defined (and variable) value of the heat transfer coefficient. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]
L.I. Rubinstein, The Stefan Problem, AMS, Providence, 1971. A.M. Meirmanov, The Stefan Problem, Walter de Gruyter, Berlin, 1992. N.L. Goldman, Inverse Stefan Problem, Kluwer, Dordrecht, 1997. N. Al-Khalidy, Application of optimization methods for solving inverse phase-change problems, Numer. Heat Transf., B Fundam. 31 (1997) 477–497. D.D. Ang, A.P.N. Dinh, D.N. Thanh, Regularization of an inverse two-phase Stefan problem, Nonlinear Anal. 34 (1998) 719–731. R. Grzymkowski, D. Słota, One-phase inverse Stefan problems solved by Adomian decomposition method, Comput. Math. Appl. 51 (2006) 33–40. P. Jochum, The numerical solution of the inverse Stefan problem, Numer. Math. 34 (1980) 411–429. J. Liu, B. Guerrier, A comparative study of domain embedding methods for regularized solutions of inverse Stefan problems, Int. J. Numer. Methods Eng. 40 (1997) 3579–3600. D. Słota, Influence of the mutation operator on the solution of an inverse Stefan problem by genetic algorithms, in: V.N. Alexandrov, et al., (Eds.), Computational Science — ICCS 2006, LNCS 3991, Springer-Verlag, Berlin, 2006, pp. 786–789. D. Słota, Solving the inverse Stefan design problem using genetic algorithms, Inverse Probl. Sci. Eng. (in review). V.R. Voller, Enthalpy method for inverse Stefan problems, Numer. Heat Transf., B Fundam. 21 (1992) 41–55. N. Zabaras, Y. Ruan, O. Richmond, Design of two-dimensional Stefan processes with desired freezing front motions, Numer. Heat Transf., B Fundam. 21 (1992) 307–325. N. Zabaras, S. Kang, On the solution of an ill-posed design solidification problem using minimization techniques in finite- and infinitedimensional function space, Int. J. Numer. Methods Eng. 36 (1993) 3973–3990. N. Zabaras, K. Yuan, Dynamic programming approach to the inverse Stefan design problem, Numer. Heat Transf., B Fundam. 26 (1994) 97–104. M. Slodička, H.D. Schepper, Determination of the heat-transfer coefficient during soldification of alloys, Comput. Methods Appl. Mech. Eng. 194 (2005) 491–498. R. Grzymkowski, D. Słota, Determination of unknown thermal coefficient in continuous casting, in: H.G. Fritz (Ed.), Proc. of the 3rd ESAFORM Conference on Material Forming, Universität Stuttgart, Stuttgart, 2000, pp. X11–X14. R. Grzymkowski, D. Słota, Numerical calculations of the heat-transfer coefficient during solidification of alloys, in: B. Sarler, C.A. Brebbia (Eds.), Moving Boundaries VI, Wit Press, Southampton, 2001, pp. 41–50. R. Grzymkowski, D. Słota, Numerical method for multi-phase inverse Stefan design problems, Arch. Metall. Mater. 51 (2006) 161–172. R. Grzymkowski, D. Słota, Multi-phase inverse Stefan problems solved by approximation method, in: R. Wyrzykowski, et al., (Eds.), Parallel Processing and Applied Mathematic, LNCS 2328, Springer-Verlag, Berlin, 2002, pp. 679–686. K. Kurpisz, A.J. Nowak, Inverse Thermal Problems, Comput. Mech. Publ., Southampton, 1995. A.N. Tikhonov, V.Y. Arsenin, Solution of Ill-Posed Problems, Wiley & Sons, New York, 1977. D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Reading, 1989. Z. Michalewicz, Genetic Algorithms + Data Structures = Evolution Programs, Springer-Verlag, Berlin, 1996. A. Osyczka, Evolutionary Algorithms for Single and Multicriteria Design Optimization, Physica-Verlag, Heidelberg, 2002. A. Kapusta, B. Mochnacki, The analysis of heat transfer processes in the cylindrical radial continuous casting volume, Bull. Pol. Acad. Sci., Tech. Sci. 36 (1988) 309–320. E. Majchrzak, B. Mochnacki, Application of the BEM in the thermal theory of foundry, Eng. Anal. Bound. Elem. 16 (1995) 99–121. J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1996. B. Furenes, B. Lie, Using event location in finite-difference methods for phase-change problems, Numer. Heat Transf., B Fundam. 50 (2006) 143–155. Rizwan-uddin, One-dimensional phase change with periodic boundary conditions, Numer. Heat Transf., A Appl. 35 (1999) 361–372. M. Zerroukat, C. R. Chatwin, Computational Moving Boundary Problems, Research Studies Press, Taunton, 1994. J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313.