- Email: [email protected]
Restoring boundary conditions in the solidification of pure metals
Computers and Structures 89 (2011) 48–54
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Restoring boundary conditions in the solidification of pure metals q Damian Słota ⇑ Institute of Mathematics, Silesian University of Technology, Kaszubska 23, Gliwice 44-100, Poland
a r t i c l e
i n f o
Article history: Received 3 December 2009 Accepted 5 August 2010 Available online 1 September 2010 Keywords: Solidification Inverse Stefan problem Genetic algorithms Continuous casting
a b s t r a c t The scope of the paper is an algorithm reconstructing the boundary conditions (heat flux and heat transfer coefficient) in the solidification of pure metals on the grounds of temperature measurements. For the verification of the algorithm experimental data derived in the course of the solidification of aluminum were used. An example of the application of the algorithm for designating the cooling conditions in continuous casting is provided, when the values of temperature at selected points on the boundary of the casting are known. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The distribution of temperature in objects exchanging heat with the surroundings may be designated from the equation or a system of heat transfer equations with all the parameters appearing in the equations and on the grounds of the boundary conditions and the initial condition for unsteady states. The task formulated in such way is a direct problem. In other words, direct problems involve the designation of the effects, i.e. temperature distribution, on the bases of the stimuli caused by the phenomenon of heat convection, the form of which is determined by the boundary conditions. Tests run on real objects as well as design tasks often enforce the inversion of the problem, for example, when the thermal responses of the system are known, in form of, for example, temperature fields on specific surfaces inside the object, temperature at measurement points, but when all the parameters evoking such responses are not completely given. In such tasks, the unknowns, apart from the function describing the temperature distribution, are other functions or parameters: functions describing the boundary conditions, thermophysical parameters, geometrical dimensions of the objects, and others. Such tasks are referred to as inverse problems and constitute very useful tools for the analyses of heat exchange processes, including the solidification of casts [1–3]. An essential case of technical problems involved in the heat transfer processes is Stefan problem [4]. This problem entails mathematical models describing thermal processes that are char-
q Extended form of paper presented at the ECCOMAS International Seminar on Inverse Methods of Mechanics, IPM 2009, April 23–25, 2009, Rzeszów-Lancut, Poland. ⇑ Tel./fax: +48 032 2372864. E-mail address: [email protected]
0045-7949/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2010.08.002
acterized by phase changes. They include: solidification of pure metals, growth of crystals, freezing of food, freezing of the ground, etc. The Stefan problem involves a simultaneous determination of the temperature distribution and the location of the moving boundary (freezing front) separating the sub-domains occupied by the liquid and solid phase. Whereas, in the inverse Stefan problem it is often assumed that additional information compensating the absence of input data is the knowledge of the location of the moving boundary, its velocity towards the normal direction or the temperature at selected points of the domain. Tasks utilising additional information on the location of the moving boundary are often referred as design problems. The inverse Stefan problem was discussed in detail in numerous works published by Zabaras and his colleagues. In [5,6] the considered case entails additional information on the temperature at selected points of the solid phase, whereas the location and velocity of the boundary of the phase change [5] or the heat flux at the boundary [6] must be found. The solution takes advantage of the deforming finite elements method, which assumes that the nodes of the finite elements change their location together with the motion of the phase change front. The minimization of the error in view of the calculated and given temperature of the phase change front is used. In [7] the inverse problem was tested in the domain where in the liquid phase heat exchange occurs by means of conduction and convection. This task was divided into two independent problems, the first one of which is a direct convection problem in the liquid phase, whereas the second one is the inverse problem in the solid phase. Yet, in [8] both component problems are treated as inverse. Similar issues were also considered in [9]. In [10,11] the authors used the Adomian decomposition method and variational iteration method in conjunction with the optimization to the approximate solution of one-phase inverse Stefan problem. Whereas in [12,13] the inverse design Stefan problem was
49
D. Słota / Computers and Structures 89 (2011) 48–54
considered, where heat transfer coefficient was designated on the boundary of the domain. A numerical solution of two-phase Stefan problem formulated for continuous casting process was given in [14,15]. In [16] however, the authors designated the location of the phase change front in continuous casting by solving an inverse geometry problem. Nawrat and Skorek in [17] used the Stefan problem for designating the thermal resistance of the gap between the casting and the crystallizer during the process of continuous casting. In [18,19] the inverse design problem involves the determination of a cooling strategy for the secondary cooling zone (the water sprays) in continuous casting process. The method of determining the optimal cooling strategy for the secondary cooling zone in continuous casting process with changing casting speed is presented in paper [20]. The evaluation of the heat transfer coefficients along the secondary cooling zone in the continuous casting of steel is presented in paper [3], while in [21,22] the selection of the cooling strategy is for the primary and secondary cooling zone. In this paper, the reconstruction of the function describing the boundary conditions (heat flux and heat transfer coefficient) in the solidification of pure metals is presented. An algorithm will be presented that enables the solution of the inverse Stefan problem, where additional information consists of temperature measurements at selected points of the domain. In numerical calculations a genetic algorithm [23] and alternating phase truncation methods were used [13]. The algorithm was implemented in C. To verify the algorithm, experimental data obtained in the course of the solidification of aluminum were used. Also, an example of using the derived algorithm for designating the cooling conditions in continuous casting is given. In this example the temperature measurements on the boundary of the casting (obtained, for example, from the infrared camera) are used. 2. Problem formulation The Stefan problem is a mathematical model of the solidification of pure metals [4]. It involves the designation of temperature distribution and the location of the phase changes boundary, when the initial condition, boundary conditions and physical parameters of the material are known. Let us assume that the solidifying metal is contained in domain X Rn , the boundary of which is surface C = oX. Domain X is a sum of two sub-domains X1(t) and X2(t) that change in time (Fig. 1). Let X1(t) denote this sub-domain which is occupied by the liquid phase and X2(t) the sub-domain occupied by the solid phase. Let us assume that the changing in time smooth surface Cg(t) (in shorter form also expressed as Cg) separates domains X1(t) and X2(t). We shall consider a case, in which the temperature distribution Tk, in each sub-domain Xk(t), fulfills the heat transfer equation [4]:
c k .k
@T k ðr; tÞ ¼ divðkk rT k ðr; tÞÞ @t
ð1Þ
for r 2 Xk(t), k = 1, 2, t 2 [0, t*), where ck, .k and kk are the specific heat, the mass density and the thermal conductivity in the liquid phase (k = 1) and solid phase (k = 2), respectively. Let us assume that at the initial moment (t = 0) the whole domain is occupied by the liquid phase. Therefore the initial condition has the form (T0 > T*):
T 1 ðr; 0Þ ¼ T 0
ð2Þ
for r 2 X1(0) = X. On the boundary C = C1(t) [ C2(t) [ C3(t) (hereunder, referred to as: C1, C2 and C3) of domain X the boundary conditions are given. On the boundary C1 a homogeneous boundary condition of the second kind is designated:
@T k ðr; tÞ ¼ 0 @n
ð3Þ
for r 2 C1(t). On the boundary C2 [ C3 the boundary condition of the second and third kind are designated:
kk
@T k ðr; tÞ ¼ qðr; tÞ @n
ð4Þ
for r 2 C2, and
kk
@T k ðr; tÞ ¼ aðr; tÞðT k ðr; tÞ T 1 Þ @n
ð5Þ
for r 2 C3, where q is the heat flux, a is the heat transfer coefficient and T1 is the surrounding temperature. Whereas, on the freezing front Cg the condition of temperature continuity and the Stefan condition are designated (r 2 Cg):
T 1 ðrÞ ¼ T 2 ðrÞ ¼ T ;
ð6Þ
L.2 v n
ð7Þ
@T 1 ðrÞ @T 2 ðrÞ ¼ k1 þ k ; 2 @n Cg @n Cg
where T* is the temperature of the phase change, L is latent heat of solidification, v n is the freezing front velocity vector in a normal direction. The discussed inverse Stefan problem consists in finding functions q and a describing the heat flux and the heat transfer coefficient on the boundaries C2 and C3, respectively. As some additional information in the inverse Stefan problem we have the temperature measurements:
T 2 ðr i ; t j Þ ¼ U ij
ð8Þ
for i = 1, 2, . . ., N1 and j = 1, 2, . . ., N2, where N1 denotes the number of the sensors, N2 denotes the number of the measurements from each sensor. Accordingly, we designate function f(r, t), describing the boundary conditions, in the following way:
f ðr; tÞ ¼
qðr; tÞ;
for r 2 C2 ;
aðr; tÞ; for r 2 C3 :
ð9Þ
Using the given temperatures Uij and calculated temperatures Tij (for a given form of functions q and a), we may construct the following functional:
Jðf Þ ¼
N1 X N2 X ðT ij U ij Þ2 : i¼1
ð10Þ
j¼1
The best solution of our problem corresponds to the minimum of this functional [24]. To find the minimum of the functional (10) a genetic algorithm is used. 3. Genetic algorithm
Fig. 1. Schematic domain of the problem.
The genetic algorithm is a method of optimization inspired by the biological mechanisms of evolution [23]. Genetic algorithms transform a set of chains with settled length that are called chromosomes (representing the vector of decision variables). A single element of the chain is called a gene. In a classic genetic algorithm
50
D. Słota / Computers and Structures 89 (2011) 48–54
the genes assume one of the two values: either zero or one, which is the so called binary coding. Also, real coding (real number representation) is used, where the chromosomes are chains with as many components as variables of the objective function, and each gene takes real values from an appropriate range. It is such type of coding that was used in this paper. The genetic algorithm starts with a population of individuals (chromosomes) randomly generated in a defined search space. Next, the value of the fitting function is calculated for each chromosome, describing whether the vector of decision variables represented by a given chromosome is a good approximation of the sought solution. The fitting function my be equal to the objective function (10) (which fact is used by the author) or dependent on it in another manner [23]. On the grounds of the fitting function values, the chromosomes are selected for further processing. The selection consists in choosing a set of chromosomes from the current population. The set that is selected has the same size as the initial population (some chromosomes may be selected many times). The tournament selection was applied in the calculation. This selection is carried out in such manner that two chromosomes are drawn and the one with better fitness goes to a new generation. There are as many draws as the individuals that the new generation is supposed to include. Because the selection cannot introduce new individuals, i.e. new points in the search space, two other mechanisms are used to achieve this: crossover and mutation. Crossover combines the genomes of the selected parents to produce an offspring. In this paper arithmetical crossover was applied as the crossover operator, where as a result of crossing two chromosomes x1 ¼ x11 ; x12 ; . . . ; x1n and x2 ¼ x21 ; x22 ; . . . ; x2n , their linear combinations are obtained: 0
x1 ¼ rx1 þ ð1 rÞx2 ; 0
x2 ¼ rx2 þ ð1 rÞx1 ;
xi þ Dðs; xui xi Þ; xi Dðs; xi xli Þ
HðTÞ ¼
Z
T
cðuÞ.ðuÞdu þ gðTÞL.2 ;
0
ð15Þ
where
gðTÞ ¼
1 for T > T ; 0
for T 6 T :
ð16Þ
Function H(T) is discontinuous at the point given by the temperature of the phase change T*. Its left-hand and right-hand limits at this point will be denoted as Hs and Hl:
Hs ¼
Z
T
cðuÞ.ðuÞdu;
ð17Þ
0
Hl ¼ Hs þ L.2 :
ð18Þ
ð11Þ ð12Þ
V 1 ðr; t i Þ ¼ maxfHl ; Hðr; t i Þg:
ð13Þ
and a decision is taken at random which from the above formulas should be applied, where:
s Dðs; xÞ ¼ x 1 r ð1NÞbm
To solve a direct Stefan problem, the alternating phase truncation method was applied [13]. In this method in place of temperature T we insert an enthalpy:
If we use Eq. (15) in the Stefan problem, we will obtain, in both phases, the heat conduction equation where the temperature will be replaced with enthalpy. The algorithm of the alternating phase truncation method (for one time’s step) consists of two stages. Let’s assume that we know the distribution of enthalpy H(r, ti) in time ti from the initial condition or from the previous step of the calculations. In the first stage, we reduce the entire domain to the liquid phase, i.e. to the points at which the value of the enthalpy is smaller than Hl, we supply (conventionally) such quantity of heat that the enthalpy equals to Hl. The pseudo-initial condition for the first stage of the calculations has the form:
where parameter r is a random number with a uniform distribution from the domain [0, 1]. The mutation operator modifies, with given probability, a random genome from a random chromosome. In the calculations, a non-uniform mutation operator was used as well. During mutation, the xi gene from chromosome x = (x1, . . . , xi, . . . , xn) is transformed according to the equation:
x0i ¼
4. Alternating phase truncation method
ð14Þ
and r is a random number with a uniform distribution from the domain [0, 1], s is the current generation number, N is the maximum number of generations and bm is a constant parameter (in the calculations, bm = 2 was assumed). In this paper the elitist model was also applied in the algorithm. In the elitist model 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. The genetic algorithm is iterated, and a new population with constant size generated, as long as the value of the fitness function of the best individual in the population does not change significantly over several generations, or, as long as the maximal number of generations is obtained.
The so obtained heat transfer problem in a one-phase domain can be solved by one of the known methods (e.g. the finite difference method or the finite elements method), thereby obtaining an b 1 ðr; tiþ1 Þ. At points to which approximate distribution of enthalpy: V we have supplied a certain amount of heat, the same amount must be now deducted:
b 1 ðr; tiþ1 Þ þ ðHðr; t i Þ V 1 ðr; t i ÞÞ: V 1 ðr; t iþ1 Þ ¼ V After this operation we obtain the distribution of enthalpy V1(r, ti+1), which is treated as a starting point for the second stage of calculations (at moment ti). In the second stage, we reduce the whole domain to the solid phase, i.e. at those points of the domain where the enthalpy value is higher than Hs, we take out (conventionally) such amount of heat that would allow the enthalpy to adopt the value equal to Hs. Thus, the pseudo-initial condition for this stage of calculations takes the following form:
V 2 ðr; t i Þ ¼ minfHs ; V 1 ðr; t iþ1 Þg: Like in the first stage, we find an approximate distribution of enb 2 ðr; t iþ1 Þ. At the end of the second stage, at the points where thalpy V we artificially take out a certain amount of heat, we add the same amount of heat:
b 2 ðr; t iþ1 Þ þ ð V b 1 ðr; t iþ1 Þ V 2 ðr; ti ÞÞ: Hðr; tiþ1 Þ ¼ V The computation of the distribution of enthalpy H(r, ti+1) completes the second stage and at the same time, one step of the calculations (transfer from time ti to time ti+1) of the alternating phase truncation method. In the alternating phase truncation method, for each time step, the heat transfer equation is solved twice. Therefore, we must take
51
D. Słota / Computers and Structures 89 (2011) 48–54
into appropriate consideration the boundary conditions, so that they would influence the discussed system only over time Dt, and not 2 Dt. In the first stage of the alternating phase truncation method, we take account of the real boundary conditions only on those boundary fragments where the liquid phase contacts the surroundings. At the same time, we isolate the remaining fragments of the boundary. In the second stage, we take into account the real boundary conditions on those boundary fragments where the solid phase contacts the surroundings and we isolate the remaining fragments of the boundary. 5. Experimental verification To verify the derived algorithm experimental data obtained in the solidification of aluminum EN AW-Al99.5 were used. The experiment was conducted by means of the UMSA (Universal Metallurgical Simulator and Analyzer) designed for analysing thermal processes occurring in metals [25]. Four cylinder samples with the diameter of 18 [mm] and height of 20 [mm] were used. The material was melted in an induction crucible furnace and cast into a 25 [mm] diameter graphite mould. Next, it was mechanically treated to the required dimensions. In the first two samples the thermocouple was located on the sample axis, whereas in the next two at the distance of 4.5 [mm] from the sample axis. The bottom and top surfaces of the samples were thermally insulated. Three full loops of the melting and solidification of the material were carried out, i.e. three temperature distributions were obtained from each sample. In view of the geometry of the domain and thermal symmetry, the modeling involved a two-phase, axisymmetrical one dimensional Stefan problem. The calculations reconstructed only function a (thus, C2 = ; was assumed), describing the heat transfer coefficient, dependent on a different number of parameters:
aðtÞ ¼ aðt; a1 ; . . . ; an Þ; n 2 f1; 3; 6; 10; 15g: To approximate the heat transfer coefficient Bezier’s curves were used [26] (in the case of one parameter, the heat transfer coefficient was approximated by a constant function). The following values of the parameters describing the process were assumed: b = 0.009 [m], k1 = 104 [W/(mK)], k2 = 240 [W/(mK)], c1 = 1290 [J/(kg K)], c2 = 1000 [J/(kg K)], .1 = 2380 [kg/m3], .2 = 2679 [kg/m3], L = 390,000 [J/kg], solidification temperature T* = 930 [K], surrounding temperature T1 = 298 [K] and initial temperature T0 = 1013 [K]. Moreover, the following values of the parameters of the genetic algorithm were used: population size npop = 70, crossover probability pc = 0.7 and mutation probability pm = 0.1. The number of generations N was changed together with increasing the number of
Fig. 2. The heat transfer coefficient reconstructed for a different number of parameter.
the sought parameters, so, accordingly, for one and three parameters the number of generations was 100, for six and ten parameters it was 1000, whereas for fifteen it equalled 2000. The operators and the values of the genetic algorithm parameters applied in the calculations were selected on the basis of a number of numerical experiments carried out for an design inverse Stefan problem [27–29]. In the alternating phase truncation method the finite differences method was used and the calculation run on the grid with discrete steps equal to Dt = 0.1 and Dx = b/500. In Fig. 2 the heat transfer coefficient reconstructed for different numbers of the sought parameters is shown, in the case where the thermocouple was located on the axis of the sample. In Table 1 mean and maximal relative and absolute errors of reconstructing the cooling curve are compiled for the heat transfer coefficients derived for different number of the parameters. Whereas in Fig. 3 the cooling curve designated in the course of the experiment is shown as well as its reconstruction for the heat transfer coefficient designated for fifteen parameters. Absolute errors of reconstructing the cooling curve are shown in Fig. 4. The next two figures illustrate the impact of the maximal number of generations N (used in the genetic algorithm) on the shape of the curve describing the heat transfer coefficient, and errors in recreating the cooling curve. Fig. 5 shows changes in the shape of the curve describing the heat transfer coefficient depending on the maximal number of generations (N = 200, 500, 1000, 2000). Fig. 6 shows the impact of the maximal number of generations on the size of errors in reconstructing the measured cooling curve. In this case, mean relative percentage errors decrease from 0.33% for N = 200 to 0.12% for N = 2000. But, mean absolute errors decrease from 2.5 [K] to 0.91 [K]. Likewise, the maximal values of errors, decrease from 2.89% to 1.41% and from 27.79 [K] to 12.94 [K], respectively. As indicated by the above mentioned results, an increase in the number of parameters leads to better reconstruction of the solidification curve. At the same time, the mean and maximal values of both relative and absolute errors decrease. Further increase in the number of parameters does not; however, result in considerable Table 1 Errors in the reconstruction of the cooling curve (dm – mean relative percentage error, dmax – maximal relative percentage error, Dm – mean absolute error, Dmax – maximal absolute error). n
1
3
6
10
15
dm [%] dmax [%] Dm [K] Dmax [K]
4.84 10.54 27.84 68.91
0.56 4.21 3.94 40.97
0.34 1.91 2.44 17.46
0.23 1.47 1.52 13.51
0.12 1.41 0.91 12.94
Fig. 3. The cooling curve (solid line – measurement data, dots – the curve reconstructed for fifteen parameters).
52
D. Słota / Computers and Structures 89 (2011) 48–54
We shall discuss a 3-D model of a vertical machine for continuous casting, working in an undisturbed cycle. Let us further assume that the heat transfer takes place only in the direction perpendicular to the casting’s axis. Such assumption stems from the fact that the quantity of heat transferred towards the direction of the casting motion is negligible in comparison with the quantity transferred towards the direction perpendicular to the casting’s axis [30]. In view of the above assumptions and the thermal symmetry, the heat exchange process, including pseudo-steady temperature field and the position of the phase change boundary, is described by a two-phase Stefan problem, where time t and casting velocity w are bound by the following relation: t = z/w, where z is a spatial variable along the casting axis. The boundary of the casting domain was divided into seven parts (Fig. 7): Fig. 4. The absolute errors of the reconstruction of the cooling curve for fifteen parameters.
C0 ¼ fðx; y; 0Þ; x 2 ½0; b; y 2 ½0; dg; C11 ¼ fð0; y; zÞ; y 2 ½0; d; z 2 ½0; z g; C12 ¼ fðx; 0; zÞ; x 2 ½0; b; z 2 ½0; z g; C21 ¼ fðb; y; zÞ; y 2 ½0; d; z 2 ½0; z1 g; C22 ¼ fðx; d; zÞ; x 2 ½0; b; z 2 ½0; z1 g; C31 ¼ fðb; y; zÞ; y 2 ½0; d; z 2 ðz1 ; z g; C32 ¼ fðx; d; zÞ; x 2 ½0; b; z 2 ðz1 ; z g; where the initial and boundary conditions will be designated. We seek the heat flux in the continuous cast mould (C2 = C21 [ C22) and the heat transfer coefficient in the secondary cooling zone (C3 = C31 [ C32). Now, we designate function f(r), describing the boundary conditions, in the following way:
f ðrÞ ¼
Fig. 5. The heat transfer coefficient reconstructed for fifteen parameters and different number of generations.
qj ðrÞ;
for r 2 C2j ;
j 2 f1; 2g;
aj ðrÞ; for r 2 C3j ; j 2 f1; 2g:
We shall now present the case where a three-dimensional model of continuous casting of copper is described by the following parameter [14]: b = 0.1 [m], d = 0.08 [m], z1 = 0.22 [m], k1 = k2 = 370 [W/(m K)], c1 = c2 = 400 [J/(kg K)], .1 = .2 = 8900 [kg/m3], L = 200,000 [J/kg], w = 0.002 [m/s], T* = 1356 [K], T1 = 323 [K] and T0 = 1423 [K]. The exact cooling conditions were described by functions: q1(r) = 350,000, q2(r) = 400,000, a1(r) = 3500, a2(r) = 4000. It was
Fig. 6. Errors in the reconstruction of the cooling curve for fifteen parameters and different number of generations (dm – mean relative percentage error, Dm – mean absolute error).
improvement of the reconstruction, yet, prolonged the calculation time. Thus, the increase in the number of the sought parameters requires increased maximal number of generations. As far as other measurements are concerned, similar conclusions may be drawn. 6. Continuous casting Let us now give an example of applying the discussed algorithm for designating the cooling conditions in continuous casting, when the values of temperature at selected points on the boundary of the casting are known (obtained from the infrared camera).
Fig. 7. Domain of the continuous casting process.
D. Słota / Computers and Structures 89 (2011) 48–54
assumed in the calculations that the values of temperature are known at one, two, three or four points located at the boundary of the domain. Two points are located on the boundary for which y = d at the distance of 7 [mm] and 17 [mm] from the point (b, d) (these points shall be marked as C and A, respectively). Whereas the remaining two points are located on the boundary for which x = b at the distance of 8 [mm] and 18 [mm] from the point (b, d) (they shall be marked as D and B, respectively). The exacted values of temperature and the values burdened with random error with normal distribution equalling 1% and 2% were used in the calculations. The following values of the genetic algorithm were used in the calculations: population size npop = 70, maximal number of generations N = 500, crossover probability pc = 0.7 and mutation probability pm = 0.1. In Table 2 the results of reconstructing the sought boundary conditions are compiled, together with relative percentage errors of such reconstruction. The results were derived for variable number of measurement points (N1, from four ABCD to one A) and exact initial data. Clearly, each time the boundary conditions are reconstructed with minimal errors, which is a consequence of the assumed criterion of finishing the algorithm. For exact input data
Table 2 Results of the reconstruction of the boundary conditions for different number of measurement points. f
Error [%]
f
Error [%]
ABCD 350014.23 399985.60 3499.62 4000.23
0.004 0.004 0.011 0.007
ABC 350146.01 399848.58 3498.22 4001.92
0.042 0.038 0.051 0.048
AB 350004.24 399972.12 3499.33 4000.82
0.001 0.007 0.019 0.020
A 349906.98 400081.57 3502.54 3996.59
0.027 0.020 0.073 0.085
53
we receive the solutions for which the value of the objective function do not exceed 0.03. In Fig. 8 the errors of reconstructing the sought heat fluxes and heat transfer coefficients concern the case when the input data are burdened with errors. It is evident that each time the errors in the reconstruction of the boundary conditions (with distorted data) are considerably smaller than the errors in the input data. For the disturbance of 2% the errors do not exceed 1.07%, whereas for 1% the errors are below 0.25%. The figure renders the results for three measurement points (ABC), in other cases, similar results were obtained.
Fig. 9. The absolute errors of the reconstruction of the temperature at the measurement point B (calculations for three measurement points ABC and exact input data).
Fig. 10. The absolute errors of the reconstruction of the temperature at the measurement point A (calculations for three measurement points ABC and distorted input data by 2%).
Fig. 8. Errors in the reconstruction of the cooling conditions for three measurement points (ABC) and different size of the errors in the input data.
Table 3 Errors in the reconstruction of the cooling curve in measurement points ABC (dm – mean relative percentage error, dmax – maximal relative percentage error, Dm – mean absolute error, Dmax – maximal absolute error). Per.
0%
1%
2%
dm [%] dmax [%] Dm [K] Dmax [K]
0.003 0.010 0.004 0.025
0.017 0.050 0.034 0.170
0.080 0.270 0.138 0.883
Fig. 11. Errors in the reconstruction of the cooling conditions for three measurement points (ABC) and different number of temperature measurements (25, 50 and 100).
54
D. Słota / Computers and Structures 89 (2011) 48–54
The temperatures at the measuring points were reconstructed very well each time, with maximal errors not exceeding 0.01% (0.025 [K]) for exact input data and 0.27% (0.89 [K]) for data burdened by errors of 2% (see Table 3). Figs. 9 and 10 present the absolute errors in the reconstruction of temperature at the measurement point for exact and distorted input data. At the same time, in Fig. 11 errors in the reconstruction of the heat fluxes and heat transfer coefficients were compiled for different number of temperature measurements (N2 2 {25, 50, 100}). As may be observed, a decrease in the number of measurements does not lead to a significant increase in the error of reconstructing the boundary conditions, the maximal error for 50 measurements was 0.08%, whereas for 25 measurements 0.11%. 7. Conclusion The scope of the paper is the algorithm for reconstructing the boundary conditions in an inverse Stefan problem. To verify the algorithm experimental data obtained in the process of the solidification of aluminum were used. In the experiment twelve cooling curves were derived. An attempt at modeling the curves by means of two-phase axisymmetrical one dimensional Stefan problem rendered very good results, which prove that even such simple mathematical model of metal solidification in conjunction with the algorithm provided very good reconstruction of the experimental data. Exemplary application of the algorithm was also shown for designating the cooling conditions of continuous casting when the values of temperature at selected points on the boundary of the casting are known (by means of an infrared camera). The featured examples of calculations show a very good approximation of the exact solution and stability of the algorithm in terms of the input data errors. For exact input data, the boundary conditions are reconstructed, in each case, with minimal errors. Whereas, for distorted input data, the errors in the reconstructed boundary conditions are, in each case, smaller than the errors in the input data. Furthermore, smaller number of temperature measurements does not result in significant changes in the reconstructed boundary conditions. Comparison of the results obtained by using the genetic algorithms and the Nelder–Mead method for design inverse Stefan problem was presented in papers [12,13]. The application of genetic algorithms enhanced the accuracy of the results obtained. Continuous casting is nowadays a widely used method of the production of casts. A proper design of the cooling system exerts an influence on securing suitable quality of the produced casts. Thus, appropriate selection of a cooling strategy for the cooling zones is of big importance. The presented method can be easy applied to solve design problems of different types, e.g. for the design of continuous casting installations (incl. the selection of the length of secondary cooling zones, the number of jets installed in individual zones, etc.). Acknowledgements I wish to thank the reviewers for their valuable criticisms and suggestions, leading to the present improved version of my paper. I would like to express my gratitude to Mirosława Pawlyta, Waldemar Kwas´ny and Mariusz Król for their assistance in obtaining the experimental data. The research was financed by the resources allocated in years 2007–2009 under research project No. N N512 3348 33. The calculations were made at the Warsaw University
Interdisciplinary Centre for Mathematical and Computational Modelling (ICM) under computational Grant No. G30-7. References [1] Dvorkin EN, Cavaliere MA, Goldschmit MB. Finite element models in the steel industry. Part I: Simulation of flat product manufacturing processes. Comput Struct 2003;81:559–73. [2] Imani A, Ranjbar AA, Esmkhani M. Simultaneous estimation of temperaturedependent thermal conductivity and heat capacity based on modified genetic algorithm. Inverse Probl Sci Eng 2006;14:767–83. [3] Santos CA, Garcia A, Frick CR, Spim JA. Evaluation of heat transfer coefficient along the secondary cooling zones in the continuous casting of steel billets. Inverse Probl Sci Eng 2006;14:687–700. [4] Gupta SC. The classical Stefan problem. Basic concepts, modelling and analysis. Amsterdam: Elsevier; 2003. [5] Zabaras N, Ruan Y. A deforming finite element method analysis of inverse Stefan problem. Int J Numer Meth Eng 1989;28:295–313. [6] Zabaras N. Inverse finite element techniques for the analysis of solidification processes. Int J Numer Meth Eng 1990;29:1569–87. [7] Zabaras N, Nguyen TH. Control of the freezing interface morphology in solidification processes in the presence of natural convection. Int J Numer Meth Eng 1995;38:1555–78. [8] Yang GZ, Zabaras N. An adjoint method for the inverse design of solidification processes with natural convection. Int J Numer Meth Eng 1998;42:1121–44. [9] Okamoto K, Li BQ. A regularization method for the inverse design of solidification processes with natural convection. Int J Heat Mass Transfer 2007;50:4409–23. [10] Grzymkowski R, Słota D. One-phase inverse Stefan problems solved by Adomian decomposition method. Comput Math Appl 2006;51:33–40. [11] Słota D. Direct and inverse one-phase Stefan problem solved by variational iteration method. Comput Math Appl 2007;54:1139–46. [12] Słota D. Using genetic algorithms for the determination of an heat transfer coefficient in three-phase inverse Stefan problem. Int Comm Heat & Mass Transfer 2008;35:149–56. [13] Słota D. Solving the inverse Stefan design problem using genetic algorithms. Inverse Probl Sci Eng 2008;16:829–46. [14] Fic A, Nowak AJ, Białecki R. Heat transfer analysis of the continuous casting process by the front tracking BEM. Eng Anal Bound Elem 2000;24:215–23. [15] Cholewa R, Nowak AJ, Białecki RA, Wrobel LC. Cubic Bezier splines for BEM heat transfer analysis of the 2-D continuous casting problems. Comput Mech 2002;28:282–90. [16] Nowak I, Nowak AJ, Wrobel LC. Identification of phase change fronts by Bezier splines and BEM. Int J Therm Sci 2002;41:492–9. [17] Nawrat A, Skorek J. Inverse finite element technique for identification of thermal resistance of gas-gab between the ingot and mould in continuous casting of metals. Inverse Probl Sci Eng 2004;12:141–55. [18] Engl HW, Langthaler T. Numerical solution of an inverse problem connected with continuous casting of steel. Z Oper Res 1985;29:B185–99. [19] Binder A, Engl HW, Vessella S. Some inverse problems for a nonlinear parabolic equation connected with continuous casting of steel: stability estimates and regularization. Numer Funct Anal Optim 1990;11:643–71. [20] Grever W, Binder A, Engl H, Mörhwald K. Optimal cooling strategies in continuous casting of steel with variable casting speed. Inverse Probl Eng 1997;2:289–300. [21] Constalas D, Kacˇur J, Keer RV. On optimal cooling strategy for variable-speed continuous casting. Int J Numer Meth Eng 2002;53:539–65. [22] Słota D. Identification of the cooling condition in 2-D and 3-D continuous casting processes. Numer Heat Transfer B 2009;55:155–76. [23] Osyczka A. Evolutionary algorithms for single and multicriteria design optimization. Heidelberg: Physica-Verlag; 2002. [24] Özisßik MN, Orlande HRB. Inverse heat transfer: fundamentals and applications. New York: Taylor & Francis; 2000. [25] Dobrzan´ski LA, Kasprzak W, Kasprzak M, Sokolowski JH. A novel approach to the design and optimization of aluminum cast component heat treatment processes using advanced UMSA physical simulations. J Achiev Mater Manuf Eng 2007;24:139–42. [26] Foley JD, van Dam A, Freiner SK, Hughes JF, Phillips RL. Computer graphics – Principles and practice. San Diego: Addison-Wesley; 1990. [27] Słota D. Influence of the mutation operator on the solution of an inverse Stefan problem by genetic algorithms. Lect Notes Comput Sci 2006;3991:786–9. [28] Słota D. Influence of choice of the crossover operator on the solution of an inverse Stefan problem by genetic algorithms. In: Cader A, Rutkowski L, Tadeusiewicz R, Zurada J, editors. Artificial intelligence and soft computing. Warszawa: EXIT; 2006. p. 217–23. [29] Słota D. Influence of the selection method on the solution of an inverse Stefan problem using genetic algorithms. In: del Pobil AP, editor. Proceedings of the 10th IASTED International Conference on Artificial Intelligence and Soft Computing, Anaheim: ACTA Press; 2006. p. 285–90. [30] Mochnacki B, Suchy JS. Numerical methods in computations of foundry processes. Cracow: PFTA; 1995.
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Restoring boundary conditions in the solidification of pure metals q Damian Słota ⇑ Institute of Mathematics, Silesian University of Technology, Kaszubska 23, Gliwice 44-100, Poland
a r t i c l e
i n f o
Article history: Received 3 December 2009 Accepted 5 August 2010 Available online 1 September 2010 Keywords: Solidification Inverse Stefan problem Genetic algorithms Continuous casting
a b s t r a c t The scope of the paper is an algorithm reconstructing the boundary conditions (heat flux and heat transfer coefficient) in the solidification of pure metals on the grounds of temperature measurements. For the verification of the algorithm experimental data derived in the course of the solidification of aluminum were used. An example of the application of the algorithm for designating the cooling conditions in continuous casting is provided, when the values of temperature at selected points on the boundary of the casting are known. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The distribution of temperature in objects exchanging heat with the surroundings may be designated from the equation or a system of heat transfer equations with all the parameters appearing in the equations and on the grounds of the boundary conditions and the initial condition for unsteady states. The task formulated in such way is a direct problem. In other words, direct problems involve the designation of the effects, i.e. temperature distribution, on the bases of the stimuli caused by the phenomenon of heat convection, the form of which is determined by the boundary conditions. Tests run on real objects as well as design tasks often enforce the inversion of the problem, for example, when the thermal responses of the system are known, in form of, for example, temperature fields on specific surfaces inside the object, temperature at measurement points, but when all the parameters evoking such responses are not completely given. In such tasks, the unknowns, apart from the function describing the temperature distribution, are other functions or parameters: functions describing the boundary conditions, thermophysical parameters, geometrical dimensions of the objects, and others. Such tasks are referred to as inverse problems and constitute very useful tools for the analyses of heat exchange processes, including the solidification of casts [1–3]. An essential case of technical problems involved in the heat transfer processes is Stefan problem [4]. This problem entails mathematical models describing thermal processes that are char-
q Extended form of paper presented at the ECCOMAS International Seminar on Inverse Methods of Mechanics, IPM 2009, April 23–25, 2009, Rzeszów-Lancut, Poland. ⇑ Tel./fax: +48 032 2372864. E-mail address: [email protected]
0045-7949/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2010.08.002
acterized by phase changes. They include: solidification of pure metals, growth of crystals, freezing of food, freezing of the ground, etc. The Stefan problem involves a simultaneous determination of the temperature distribution and the location of the moving boundary (freezing front) separating the sub-domains occupied by the liquid and solid phase. Whereas, in the inverse Stefan problem it is often assumed that additional information compensating the absence of input data is the knowledge of the location of the moving boundary, its velocity towards the normal direction or the temperature at selected points of the domain. Tasks utilising additional information on the location of the moving boundary are often referred as design problems. The inverse Stefan problem was discussed in detail in numerous works published by Zabaras and his colleagues. In [5,6] the considered case entails additional information on the temperature at selected points of the solid phase, whereas the location and velocity of the boundary of the phase change [5] or the heat flux at the boundary [6] must be found. The solution takes advantage of the deforming finite elements method, which assumes that the nodes of the finite elements change their location together with the motion of the phase change front. The minimization of the error in view of the calculated and given temperature of the phase change front is used. In [7] the inverse problem was tested in the domain where in the liquid phase heat exchange occurs by means of conduction and convection. This task was divided into two independent problems, the first one of which is a direct convection problem in the liquid phase, whereas the second one is the inverse problem in the solid phase. Yet, in [8] both component problems are treated as inverse. Similar issues were also considered in [9]. In [10,11] the authors used the Adomian decomposition method and variational iteration method in conjunction with the optimization to the approximate solution of one-phase inverse Stefan problem. Whereas in [12,13] the inverse design Stefan problem was
49
D. Słota / Computers and Structures 89 (2011) 48–54
considered, where heat transfer coefficient was designated on the boundary of the domain. A numerical solution of two-phase Stefan problem formulated for continuous casting process was given in [14,15]. In [16] however, the authors designated the location of the phase change front in continuous casting by solving an inverse geometry problem. Nawrat and Skorek in [17] used the Stefan problem for designating the thermal resistance of the gap between the casting and the crystallizer during the process of continuous casting. In [18,19] the inverse design problem involves the determination of a cooling strategy for the secondary cooling zone (the water sprays) in continuous casting process. The method of determining the optimal cooling strategy for the secondary cooling zone in continuous casting process with changing casting speed is presented in paper [20]. The evaluation of the heat transfer coefficients along the secondary cooling zone in the continuous casting of steel is presented in paper [3], while in [21,22] the selection of the cooling strategy is for the primary and secondary cooling zone. In this paper, the reconstruction of the function describing the boundary conditions (heat flux and heat transfer coefficient) in the solidification of pure metals is presented. An algorithm will be presented that enables the solution of the inverse Stefan problem, where additional information consists of temperature measurements at selected points of the domain. In numerical calculations a genetic algorithm [23] and alternating phase truncation methods were used [13]. The algorithm was implemented in C. To verify the algorithm, experimental data obtained in the course of the solidification of aluminum were used. Also, an example of using the derived algorithm for designating the cooling conditions in continuous casting is given. In this example the temperature measurements on the boundary of the casting (obtained, for example, from the infrared camera) are used. 2. Problem formulation The Stefan problem is a mathematical model of the solidification of pure metals [4]. It involves the designation of temperature distribution and the location of the phase changes boundary, when the initial condition, boundary conditions and physical parameters of the material are known. Let us assume that the solidifying metal is contained in domain X Rn , the boundary of which is surface C = oX. Domain X is a sum of two sub-domains X1(t) and X2(t) that change in time (Fig. 1). Let X1(t) denote this sub-domain which is occupied by the liquid phase and X2(t) the sub-domain occupied by the solid phase. Let us assume that the changing in time smooth surface Cg(t) (in shorter form also expressed as Cg) separates domains X1(t) and X2(t). We shall consider a case, in which the temperature distribution Tk, in each sub-domain Xk(t), fulfills the heat transfer equation [4]:
c k .k
@T k ðr; tÞ ¼ divðkk rT k ðr; tÞÞ @t
ð1Þ
for r 2 Xk(t), k = 1, 2, t 2 [0, t*), where ck, .k and kk are the specific heat, the mass density and the thermal conductivity in the liquid phase (k = 1) and solid phase (k = 2), respectively. Let us assume that at the initial moment (t = 0) the whole domain is occupied by the liquid phase. Therefore the initial condition has the form (T0 > T*):
T 1 ðr; 0Þ ¼ T 0
ð2Þ
for r 2 X1(0) = X. On the boundary C = C1(t) [ C2(t) [ C3(t) (hereunder, referred to as: C1, C2 and C3) of domain X the boundary conditions are given. On the boundary C1 a homogeneous boundary condition of the second kind is designated:
@T k ðr; tÞ ¼ 0 @n
ð3Þ
for r 2 C1(t). On the boundary C2 [ C3 the boundary condition of the second and third kind are designated:
kk
@T k ðr; tÞ ¼ qðr; tÞ @n
ð4Þ
for r 2 C2, and
kk
@T k ðr; tÞ ¼ aðr; tÞðT k ðr; tÞ T 1 Þ @n
ð5Þ
for r 2 C3, where q is the heat flux, a is the heat transfer coefficient and T1 is the surrounding temperature. Whereas, on the freezing front Cg the condition of temperature continuity and the Stefan condition are designated (r 2 Cg):
T 1 ðrÞ ¼ T 2 ðrÞ ¼ T ;
ð6Þ
L.2 v n
ð7Þ
@T 1 ðrÞ @T 2 ðrÞ ¼ k1 þ k ; 2 @n Cg @n Cg
where T* is the temperature of the phase change, L is latent heat of solidification, v n is the freezing front velocity vector in a normal direction. The discussed inverse Stefan problem consists in finding functions q and a describing the heat flux and the heat transfer coefficient on the boundaries C2 and C3, respectively. As some additional information in the inverse Stefan problem we have the temperature measurements:
T 2 ðr i ; t j Þ ¼ U ij
ð8Þ
for i = 1, 2, . . ., N1 and j = 1, 2, . . ., N2, where N1 denotes the number of the sensors, N2 denotes the number of the measurements from each sensor. Accordingly, we designate function f(r, t), describing the boundary conditions, in the following way:
f ðr; tÞ ¼
qðr; tÞ;
for r 2 C2 ;
aðr; tÞ; for r 2 C3 :
ð9Þ
Using the given temperatures Uij and calculated temperatures Tij (for a given form of functions q and a), we may construct the following functional:
Jðf Þ ¼
N1 X N2 X ðT ij U ij Þ2 : i¼1
ð10Þ
j¼1
The best solution of our problem corresponds to the minimum of this functional [24]. To find the minimum of the functional (10) a genetic algorithm is used. 3. Genetic algorithm
Fig. 1. Schematic domain of the problem.
The genetic algorithm is a method of optimization inspired by the biological mechanisms of evolution [23]. Genetic algorithms transform a set of chains with settled length that are called chromosomes (representing the vector of decision variables). A single element of the chain is called a gene. In a classic genetic algorithm
50
D. Słota / Computers and Structures 89 (2011) 48–54
the genes assume one of the two values: either zero or one, which is the so called binary coding. Also, real coding (real number representation) is used, where the chromosomes are chains with as many components as variables of the objective function, and each gene takes real values from an appropriate range. It is such type of coding that was used in this paper. The genetic algorithm starts with a population of individuals (chromosomes) randomly generated in a defined search space. Next, the value of the fitting function is calculated for each chromosome, describing whether the vector of decision variables represented by a given chromosome is a good approximation of the sought solution. The fitting function my be equal to the objective function (10) (which fact is used by the author) or dependent on it in another manner [23]. On the grounds of the fitting function values, the chromosomes are selected for further processing. The selection consists in choosing a set of chromosomes from the current population. The set that is selected has the same size as the initial population (some chromosomes may be selected many times). The tournament selection was applied in the calculation. This selection is carried out in such manner that two chromosomes are drawn and the one with better fitness goes to a new generation. There are as many draws as the individuals that the new generation is supposed to include. Because the selection cannot introduce new individuals, i.e. new points in the search space, two other mechanisms are used to achieve this: crossover and mutation. Crossover combines the genomes of the selected parents to produce an offspring. In this paper arithmetical crossover was applied as the crossover operator, where as a result of crossing two chromosomes x1 ¼ x11 ; x12 ; . . . ; x1n and x2 ¼ x21 ; x22 ; . . . ; x2n , their linear combinations are obtained: 0
x1 ¼ rx1 þ ð1 rÞx2 ; 0
x2 ¼ rx2 þ ð1 rÞx1 ;
xi þ Dðs; xui xi Þ; xi Dðs; xi xli Þ
HðTÞ ¼
Z
T
cðuÞ.ðuÞdu þ gðTÞL.2 ;
0
ð15Þ
where
gðTÞ ¼
1 for T > T ; 0
for T 6 T :
ð16Þ
Function H(T) is discontinuous at the point given by the temperature of the phase change T*. Its left-hand and right-hand limits at this point will be denoted as Hs and Hl:
Hs ¼
Z
T
cðuÞ.ðuÞdu;
ð17Þ
0
Hl ¼ Hs þ L.2 :
ð18Þ
ð11Þ ð12Þ
V 1 ðr; t i Þ ¼ maxfHl ; Hðr; t i Þg:
ð13Þ
and a decision is taken at random which from the above formulas should be applied, where:
s Dðs; xÞ ¼ x 1 r ð1NÞbm
To solve a direct Stefan problem, the alternating phase truncation method was applied [13]. In this method in place of temperature T we insert an enthalpy:
If we use Eq. (15) in the Stefan problem, we will obtain, in both phases, the heat conduction equation where the temperature will be replaced with enthalpy. The algorithm of the alternating phase truncation method (for one time’s step) consists of two stages. Let’s assume that we know the distribution of enthalpy H(r, ti) in time ti from the initial condition or from the previous step of the calculations. In the first stage, we reduce the entire domain to the liquid phase, i.e. to the points at which the value of the enthalpy is smaller than Hl, we supply (conventionally) such quantity of heat that the enthalpy equals to Hl. The pseudo-initial condition for the first stage of the calculations has the form:
where parameter r is a random number with a uniform distribution from the domain [0, 1]. The mutation operator modifies, with given probability, a random genome from a random chromosome. In the calculations, a non-uniform mutation operator was used as well. During mutation, the xi gene from chromosome x = (x1, . . . , xi, . . . , xn) is transformed according to the equation:
x0i ¼
4. Alternating phase truncation method
ð14Þ
and r is a random number with a uniform distribution from the domain [0, 1], s is the current generation number, N is the maximum number of generations and bm is a constant parameter (in the calculations, bm = 2 was assumed). In this paper the elitist model was also applied in the algorithm. In the elitist model 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. The genetic algorithm is iterated, and a new population with constant size generated, as long as the value of the fitness function of the best individual in the population does not change significantly over several generations, or, as long as the maximal number of generations is obtained.
The so obtained heat transfer problem in a one-phase domain can be solved by one of the known methods (e.g. the finite difference method or the finite elements method), thereby obtaining an b 1 ðr; tiþ1 Þ. At points to which approximate distribution of enthalpy: V we have supplied a certain amount of heat, the same amount must be now deducted:
b 1 ðr; tiþ1 Þ þ ðHðr; t i Þ V 1 ðr; t i ÞÞ: V 1 ðr; t iþ1 Þ ¼ V After this operation we obtain the distribution of enthalpy V1(r, ti+1), which is treated as a starting point for the second stage of calculations (at moment ti). In the second stage, we reduce the whole domain to the solid phase, i.e. at those points of the domain where the enthalpy value is higher than Hs, we take out (conventionally) such amount of heat that would allow the enthalpy to adopt the value equal to Hs. Thus, the pseudo-initial condition for this stage of calculations takes the following form:
V 2 ðr; t i Þ ¼ minfHs ; V 1 ðr; t iþ1 Þg: Like in the first stage, we find an approximate distribution of enb 2 ðr; t iþ1 Þ. At the end of the second stage, at the points where thalpy V we artificially take out a certain amount of heat, we add the same amount of heat:
b 2 ðr; t iþ1 Þ þ ð V b 1 ðr; t iþ1 Þ V 2 ðr; ti ÞÞ: Hðr; tiþ1 Þ ¼ V The computation of the distribution of enthalpy H(r, ti+1) completes the second stage and at the same time, one step of the calculations (transfer from time ti to time ti+1) of the alternating phase truncation method. In the alternating phase truncation method, for each time step, the heat transfer equation is solved twice. Therefore, we must take
51
D. Słota / Computers and Structures 89 (2011) 48–54
into appropriate consideration the boundary conditions, so that they would influence the discussed system only over time Dt, and not 2 Dt. In the first stage of the alternating phase truncation method, we take account of the real boundary conditions only on those boundary fragments where the liquid phase contacts the surroundings. At the same time, we isolate the remaining fragments of the boundary. In the second stage, we take into account the real boundary conditions on those boundary fragments where the solid phase contacts the surroundings and we isolate the remaining fragments of the boundary. 5. Experimental verification To verify the derived algorithm experimental data obtained in the solidification of aluminum EN AW-Al99.5 were used. The experiment was conducted by means of the UMSA (Universal Metallurgical Simulator and Analyzer) designed for analysing thermal processes occurring in metals [25]. Four cylinder samples with the diameter of 18 [mm] and height of 20 [mm] were used. The material was melted in an induction crucible furnace and cast into a 25 [mm] diameter graphite mould. Next, it was mechanically treated to the required dimensions. In the first two samples the thermocouple was located on the sample axis, whereas in the next two at the distance of 4.5 [mm] from the sample axis. The bottom and top surfaces of the samples were thermally insulated. Three full loops of the melting and solidification of the material were carried out, i.e. three temperature distributions were obtained from each sample. In view of the geometry of the domain and thermal symmetry, the modeling involved a two-phase, axisymmetrical one dimensional Stefan problem. The calculations reconstructed only function a (thus, C2 = ; was assumed), describing the heat transfer coefficient, dependent on a different number of parameters:
aðtÞ ¼ aðt; a1 ; . . . ; an Þ; n 2 f1; 3; 6; 10; 15g: To approximate the heat transfer coefficient Bezier’s curves were used [26] (in the case of one parameter, the heat transfer coefficient was approximated by a constant function). The following values of the parameters describing the process were assumed: b = 0.009 [m], k1 = 104 [W/(mK)], k2 = 240 [W/(mK)], c1 = 1290 [J/(kg K)], c2 = 1000 [J/(kg K)], .1 = 2380 [kg/m3], .2 = 2679 [kg/m3], L = 390,000 [J/kg], solidification temperature T* = 930 [K], surrounding temperature T1 = 298 [K] and initial temperature T0 = 1013 [K]. Moreover, the following values of the parameters of the genetic algorithm were used: population size npop = 70, crossover probability pc = 0.7 and mutation probability pm = 0.1. The number of generations N was changed together with increasing the number of
Fig. 2. The heat transfer coefficient reconstructed for a different number of parameter.
the sought parameters, so, accordingly, for one and three parameters the number of generations was 100, for six and ten parameters it was 1000, whereas for fifteen it equalled 2000. The operators and the values of the genetic algorithm parameters applied in the calculations were selected on the basis of a number of numerical experiments carried out for an design inverse Stefan problem [27–29]. In the alternating phase truncation method the finite differences method was used and the calculation run on the grid with discrete steps equal to Dt = 0.1 and Dx = b/500. In Fig. 2 the heat transfer coefficient reconstructed for different numbers of the sought parameters is shown, in the case where the thermocouple was located on the axis of the sample. In Table 1 mean and maximal relative and absolute errors of reconstructing the cooling curve are compiled for the heat transfer coefficients derived for different number of the parameters. Whereas in Fig. 3 the cooling curve designated in the course of the experiment is shown as well as its reconstruction for the heat transfer coefficient designated for fifteen parameters. Absolute errors of reconstructing the cooling curve are shown in Fig. 4. The next two figures illustrate the impact of the maximal number of generations N (used in the genetic algorithm) on the shape of the curve describing the heat transfer coefficient, and errors in recreating the cooling curve. Fig. 5 shows changes in the shape of the curve describing the heat transfer coefficient depending on the maximal number of generations (N = 200, 500, 1000, 2000). Fig. 6 shows the impact of the maximal number of generations on the size of errors in reconstructing the measured cooling curve. In this case, mean relative percentage errors decrease from 0.33% for N = 200 to 0.12% for N = 2000. But, mean absolute errors decrease from 2.5 [K] to 0.91 [K]. Likewise, the maximal values of errors, decrease from 2.89% to 1.41% and from 27.79 [K] to 12.94 [K], respectively. As indicated by the above mentioned results, an increase in the number of parameters leads to better reconstruction of the solidification curve. At the same time, the mean and maximal values of both relative and absolute errors decrease. Further increase in the number of parameters does not; however, result in considerable Table 1 Errors in the reconstruction of the cooling curve (dm – mean relative percentage error, dmax – maximal relative percentage error, Dm – mean absolute error, Dmax – maximal absolute error). n
1
3
6
10
15
dm [%] dmax [%] Dm [K] Dmax [K]
4.84 10.54 27.84 68.91
0.56 4.21 3.94 40.97
0.34 1.91 2.44 17.46
0.23 1.47 1.52 13.51
0.12 1.41 0.91 12.94
Fig. 3. The cooling curve (solid line – measurement data, dots – the curve reconstructed for fifteen parameters).
52
D. Słota / Computers and Structures 89 (2011) 48–54
We shall discuss a 3-D model of a vertical machine for continuous casting, working in an undisturbed cycle. Let us further assume that the heat transfer takes place only in the direction perpendicular to the casting’s axis. Such assumption stems from the fact that the quantity of heat transferred towards the direction of the casting motion is negligible in comparison with the quantity transferred towards the direction perpendicular to the casting’s axis [30]. In view of the above assumptions and the thermal symmetry, the heat exchange process, including pseudo-steady temperature field and the position of the phase change boundary, is described by a two-phase Stefan problem, where time t and casting velocity w are bound by the following relation: t = z/w, where z is a spatial variable along the casting axis. The boundary of the casting domain was divided into seven parts (Fig. 7): Fig. 4. The absolute errors of the reconstruction of the cooling curve for fifteen parameters.
C0 ¼ fðx; y; 0Þ; x 2 ½0; b; y 2 ½0; dg; C11 ¼ fð0; y; zÞ; y 2 ½0; d; z 2 ½0; z g; C12 ¼ fðx; 0; zÞ; x 2 ½0; b; z 2 ½0; z g; C21 ¼ fðb; y; zÞ; y 2 ½0; d; z 2 ½0; z1 g; C22 ¼ fðx; d; zÞ; x 2 ½0; b; z 2 ½0; z1 g; C31 ¼ fðb; y; zÞ; y 2 ½0; d; z 2 ðz1 ; z g; C32 ¼ fðx; d; zÞ; x 2 ½0; b; z 2 ðz1 ; z g; where the initial and boundary conditions will be designated. We seek the heat flux in the continuous cast mould (C2 = C21 [ C22) and the heat transfer coefficient in the secondary cooling zone (C3 = C31 [ C32). Now, we designate function f(r), describing the boundary conditions, in the following way:
f ðrÞ ¼
Fig. 5. The heat transfer coefficient reconstructed for fifteen parameters and different number of generations.
qj ðrÞ;
for r 2 C2j ;
j 2 f1; 2g;
aj ðrÞ; for r 2 C3j ; j 2 f1; 2g:
We shall now present the case where a three-dimensional model of continuous casting of copper is described by the following parameter [14]: b = 0.1 [m], d = 0.08 [m], z1 = 0.22 [m], k1 = k2 = 370 [W/(m K)], c1 = c2 = 400 [J/(kg K)], .1 = .2 = 8900 [kg/m3], L = 200,000 [J/kg], w = 0.002 [m/s], T* = 1356 [K], T1 = 323 [K] and T0 = 1423 [K]. The exact cooling conditions were described by functions: q1(r) = 350,000, q2(r) = 400,000, a1(r) = 3500, a2(r) = 4000. It was
Fig. 6. Errors in the reconstruction of the cooling curve for fifteen parameters and different number of generations (dm – mean relative percentage error, Dm – mean absolute error).
improvement of the reconstruction, yet, prolonged the calculation time. Thus, the increase in the number of the sought parameters requires increased maximal number of generations. As far as other measurements are concerned, similar conclusions may be drawn. 6. Continuous casting Let us now give an example of applying the discussed algorithm for designating the cooling conditions in continuous casting, when the values of temperature at selected points on the boundary of the casting are known (obtained from the infrared camera).
Fig. 7. Domain of the continuous casting process.
D. Słota / Computers and Structures 89 (2011) 48–54
assumed in the calculations that the values of temperature are known at one, two, three or four points located at the boundary of the domain. Two points are located on the boundary for which y = d at the distance of 7 [mm] and 17 [mm] from the point (b, d) (these points shall be marked as C and A, respectively). Whereas the remaining two points are located on the boundary for which x = b at the distance of 8 [mm] and 18 [mm] from the point (b, d) (they shall be marked as D and B, respectively). The exacted values of temperature and the values burdened with random error with normal distribution equalling 1% and 2% were used in the calculations. The following values of the genetic algorithm were used in the calculations: population size npop = 70, maximal number of generations N = 500, crossover probability pc = 0.7 and mutation probability pm = 0.1. In Table 2 the results of reconstructing the sought boundary conditions are compiled, together with relative percentage errors of such reconstruction. The results were derived for variable number of measurement points (N1, from four ABCD to one A) and exact initial data. Clearly, each time the boundary conditions are reconstructed with minimal errors, which is a consequence of the assumed criterion of finishing the algorithm. For exact input data
Table 2 Results of the reconstruction of the boundary conditions for different number of measurement points. f
Error [%]
f
Error [%]
ABCD 350014.23 399985.60 3499.62 4000.23
0.004 0.004 0.011 0.007
ABC 350146.01 399848.58 3498.22 4001.92
0.042 0.038 0.051 0.048
AB 350004.24 399972.12 3499.33 4000.82
0.001 0.007 0.019 0.020
A 349906.98 400081.57 3502.54 3996.59
0.027 0.020 0.073 0.085
53
we receive the solutions for which the value of the objective function do not exceed 0.03. In Fig. 8 the errors of reconstructing the sought heat fluxes and heat transfer coefficients concern the case when the input data are burdened with errors. It is evident that each time the errors in the reconstruction of the boundary conditions (with distorted data) are considerably smaller than the errors in the input data. For the disturbance of 2% the errors do not exceed 1.07%, whereas for 1% the errors are below 0.25%. The figure renders the results for three measurement points (ABC), in other cases, similar results were obtained.
Fig. 9. The absolute errors of the reconstruction of the temperature at the measurement point B (calculations for three measurement points ABC and exact input data).
Fig. 10. The absolute errors of the reconstruction of the temperature at the measurement point A (calculations for three measurement points ABC and distorted input data by 2%).
Fig. 8. Errors in the reconstruction of the cooling conditions for three measurement points (ABC) and different size of the errors in the input data.
Table 3 Errors in the reconstruction of the cooling curve in measurement points ABC (dm – mean relative percentage error, dmax – maximal relative percentage error, Dm – mean absolute error, Dmax – maximal absolute error). Per.
0%
1%
2%
dm [%] dmax [%] Dm [K] Dmax [K]
0.003 0.010 0.004 0.025
0.017 0.050 0.034 0.170
0.080 0.270 0.138 0.883
Fig. 11. Errors in the reconstruction of the cooling conditions for three measurement points (ABC) and different number of temperature measurements (25, 50 and 100).
54
D. Słota / Computers and Structures 89 (2011) 48–54
The temperatures at the measuring points were reconstructed very well each time, with maximal errors not exceeding 0.01% (0.025 [K]) for exact input data and 0.27% (0.89 [K]) for data burdened by errors of 2% (see Table 3). Figs. 9 and 10 present the absolute errors in the reconstruction of temperature at the measurement point for exact and distorted input data. At the same time, in Fig. 11 errors in the reconstruction of the heat fluxes and heat transfer coefficients were compiled for different number of temperature measurements (N2 2 {25, 50, 100}). As may be observed, a decrease in the number of measurements does not lead to a significant increase in the error of reconstructing the boundary conditions, the maximal error for 50 measurements was 0.08%, whereas for 25 measurements 0.11%. 7. Conclusion The scope of the paper is the algorithm for reconstructing the boundary conditions in an inverse Stefan problem. To verify the algorithm experimental data obtained in the process of the solidification of aluminum were used. In the experiment twelve cooling curves were derived. An attempt at modeling the curves by means of two-phase axisymmetrical one dimensional Stefan problem rendered very good results, which prove that even such simple mathematical model of metal solidification in conjunction with the algorithm provided very good reconstruction of the experimental data. Exemplary application of the algorithm was also shown for designating the cooling conditions of continuous casting when the values of temperature at selected points on the boundary of the casting are known (by means of an infrared camera). The featured examples of calculations show a very good approximation of the exact solution and stability of the algorithm in terms of the input data errors. For exact input data, the boundary conditions are reconstructed, in each case, with minimal errors. Whereas, for distorted input data, the errors in the reconstructed boundary conditions are, in each case, smaller than the errors in the input data. Furthermore, smaller number of temperature measurements does not result in significant changes in the reconstructed boundary conditions. Comparison of the results obtained by using the genetic algorithms and the Nelder–Mead method for design inverse Stefan problem was presented in papers [12,13]. The application of genetic algorithms enhanced the accuracy of the results obtained. Continuous casting is nowadays a widely used method of the production of casts. A proper design of the cooling system exerts an influence on securing suitable quality of the produced casts. Thus, appropriate selection of a cooling strategy for the cooling zones is of big importance. The presented method can be easy applied to solve design problems of different types, e.g. for the design of continuous casting installations (incl. the selection of the length of secondary cooling zones, the number of jets installed in individual zones, etc.). Acknowledgements I wish to thank the reviewers for their valuable criticisms and suggestions, leading to the present improved version of my paper. I would like to express my gratitude to Mirosława Pawlyta, Waldemar Kwas´ny and Mariusz Król for their assistance in obtaining the experimental data. The research was financed by the resources allocated in years 2007–2009 under research project No. N N512 3348 33. The calculations were made at the Warsaw University
Interdisciplinary Centre for Mathematical and Computational Modelling (ICM) under computational Grant No. G30-7. References [1] Dvorkin EN, Cavaliere MA, Goldschmit MB. Finite element models in the steel industry. Part I: Simulation of flat product manufacturing processes. Comput Struct 2003;81:559–73. [2] Imani A, Ranjbar AA, Esmkhani M. Simultaneous estimation of temperaturedependent thermal conductivity and heat capacity based on modified genetic algorithm. Inverse Probl Sci Eng 2006;14:767–83. [3] Santos CA, Garcia A, Frick CR, Spim JA. Evaluation of heat transfer coefficient along the secondary cooling zones in the continuous casting of steel billets. Inverse Probl Sci Eng 2006;14:687–700. [4] Gupta SC. The classical Stefan problem. Basic concepts, modelling and analysis. Amsterdam: Elsevier; 2003. [5] Zabaras N, Ruan Y. A deforming finite element method analysis of inverse Stefan problem. Int J Numer Meth Eng 1989;28:295–313. [6] Zabaras N. Inverse finite element techniques for the analysis of solidification processes. Int J Numer Meth Eng 1990;29:1569–87. [7] Zabaras N, Nguyen TH. Control of the freezing interface morphology in solidification processes in the presence of natural convection. Int J Numer Meth Eng 1995;38:1555–78. [8] Yang GZ, Zabaras N. An adjoint method for the inverse design of solidification processes with natural convection. Int J Numer Meth Eng 1998;42:1121–44. [9] Okamoto K, Li BQ. A regularization method for the inverse design of solidification processes with natural convection. Int J Heat Mass Transfer 2007;50:4409–23. [10] Grzymkowski R, Słota D. One-phase inverse Stefan problems solved by Adomian decomposition method. Comput Math Appl 2006;51:33–40. [11] Słota D. Direct and inverse one-phase Stefan problem solved by variational iteration method. Comput Math Appl 2007;54:1139–46. [12] Słota D. Using genetic algorithms for the determination of an heat transfer coefficient in three-phase inverse Stefan problem. Int Comm Heat & Mass Transfer 2008;35:149–56. [13] Słota D. Solving the inverse Stefan design problem using genetic algorithms. Inverse Probl Sci Eng 2008;16:829–46. [14] Fic A, Nowak AJ, Białecki R. Heat transfer analysis of the continuous casting process by the front tracking BEM. Eng Anal Bound Elem 2000;24:215–23. [15] Cholewa R, Nowak AJ, Białecki RA, Wrobel LC. Cubic Bezier splines for BEM heat transfer analysis of the 2-D continuous casting problems. Comput Mech 2002;28:282–90. [16] Nowak I, Nowak AJ, Wrobel LC. Identification of phase change fronts by Bezier splines and BEM. Int J Therm Sci 2002;41:492–9. [17] Nawrat A, Skorek J. Inverse finite element technique for identification of thermal resistance of gas-gab between the ingot and mould in continuous casting of metals. Inverse Probl Sci Eng 2004;12:141–55. [18] Engl HW, Langthaler T. Numerical solution of an inverse problem connected with continuous casting of steel. Z Oper Res 1985;29:B185–99. [19] Binder A, Engl HW, Vessella S. Some inverse problems for a nonlinear parabolic equation connected with continuous casting of steel: stability estimates and regularization. Numer Funct Anal Optim 1990;11:643–71. [20] Grever W, Binder A, Engl H, Mörhwald K. Optimal cooling strategies in continuous casting of steel with variable casting speed. Inverse Probl Eng 1997;2:289–300. [21] Constalas D, Kacˇur J, Keer RV. On optimal cooling strategy for variable-speed continuous casting. Int J Numer Meth Eng 2002;53:539–65. [22] Słota D. Identification of the cooling condition in 2-D and 3-D continuous casting processes. Numer Heat Transfer B 2009;55:155–76. [23] Osyczka A. Evolutionary algorithms for single and multicriteria design optimization. Heidelberg: Physica-Verlag; 2002. [24] Özisßik MN, Orlande HRB. Inverse heat transfer: fundamentals and applications. New York: Taylor & Francis; 2000. [25] Dobrzan´ski LA, Kasprzak W, Kasprzak M, Sokolowski JH. A novel approach to the design and optimization of aluminum cast component heat treatment processes using advanced UMSA physical simulations. J Achiev Mater Manuf Eng 2007;24:139–42. [26] Foley JD, van Dam A, Freiner SK, Hughes JF, Phillips RL. Computer graphics – Principles and practice. San Diego: Addison-Wesley; 1990. [27] Słota D. Influence of the mutation operator on the solution of an inverse Stefan problem by genetic algorithms. Lect Notes Comput Sci 2006;3991:786–9. [28] Słota D. Influence of choice of the crossover operator on the solution of an inverse Stefan problem by genetic algorithms. In: Cader A, Rutkowski L, Tadeusiewicz R, Zurada J, editors. Artificial intelligence and soft computing. Warszawa: EXIT; 2006. p. 217–23. [29] Słota D. Influence of the selection method on the solution of an inverse Stefan problem using genetic algorithms. In: del Pobil AP, editor. Proceedings of the 10th IASTED International Conference on Artificial Intelligence and Soft Computing, Anaheim: ACTA Press; 2006. p. 285–90. [30] Mochnacki B, Suchy JS. Numerical methods in computations of foundry processes. Cracow: PFTA; 1995.