Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 © 2017 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3805 Original Scientific Paper Received for review: 2016-06-15 Received revised form: 2016-11-21 Accepted for publication: 2016-11-29 Inverse Method for Controlling Pure Material Solidification in Spherical Geometry Mohamed Charifi* - Rabah Zegadi University of Setif, Institute of Optics and Precision Mechanics, Algeria In this study, we present the control of the solidification process of a phase-changing, pure material described in one-dimensional spherical geometry. We used an inverse global descent method in which the gradient and the adjoint equation are constructed in continuous variables of time and space. The control variable is the temperature at the fixed boundary of the solid domain. For the desired solidification front, the control was determined using information on the heat flux deduced by heat balance. The numerical resolution was based on a finite difference method in a physical domain with a moving grid related to the evolving solidification front with time. The developed numerical model was validated using an exact built solution. The numerical results of the control problem are presented for both the exact and noisy data cases. For the noisy data, a regularization method was applied. In the case of the exactdata, a rapid control determination was achieved except for time steps near the end. The random errors effects in bruited data were considerably reduced by regularization. Keywords: phase change, interface solid/liquid, inverse problem, spherical geometry Highlights • Control of pure material solidification in spherical geometry was studied for a solidification planar front. • An inverse global method has been used to determine the temperature at the fixed boundary of the solid domain as a control. • The numerical resolution was based on the finite difference method in a physical domain with a moving grid. • The algorithm has enabled the rapid control determination with excellent accuracy. 0 INTRODUCTION Material phase change represents a major challenge in various fields (metallurgy, heat storage systems, food conservation, etc.). Each year, more than a billion tons of metals are solidified around the world, mainly ferrous alloys (steel, cast iron), and aluminium. Moreover, storing thermal energy methods used in the food and pharmaceutical industries and air conditioning are directly related to the melting and solidification phenomena. In the case of material solidification, the kinetics of the state change and, more specifically, the geometry of the interface phase transition and its evolution over time determine the structure and properties of the final state. The determination by direct measurement of the mobile interface is generally very difficult to achieve. Inverse methods are most commonly used to simulate such phenomena of phase change. The control of the solidification front of a material by simulation, in metallurgy, for example, enables modification of its mechanical properties (hardness, toughness or mechanical strength). These properties depend mainly on the parameters used in the simulation (temperature, heat flux, the geometry of the solid/liquid interface and its velocity evolution). It is, therefore, necessary to control the solid/liquid interface in order to obtain particular desired properties. The melting or solidification of pure substances and eutectic alloys are characterized by a well-defined melting temperature. The solid and liquid phases are well differentiated by a sharp interface. In the case of mixtures and non-eutectic alloys, the phase change spreads out over a temperature range where the solid and liquid phases coexist. This two-phase mushy region is limited by an interface with the liquid phase at the liquidus temperature and another interface with the solid phase at the solidus temperature [1]. Solid/liquid phase change problems have been widely discussed by several authors using either inverse or direct methods. Zabaras and Ruan [2] sequentially solved a one-dimensional inverse Stefan solidification problem. They used a deformable finite element method to calculate the position and the speed of displacement of the interface and information on the temperature measured by two or more sensors located within the solid phase. Zabrasand Kang [3] treated an iterative resolution numerical simulation of a freezing front control problem in a linear case. Samaï et al. [4] identified the position of the solidification front using an iterative descent method. Jiang and al. [5] solved an inverse problem by the conjugate gradient method using finite differences to determine the historic heat flux and the final temperature distribution. Tikhonov's zero order regularization method was introduced to stabilize the inverse solution. Hetmaniok and Slota *Corr. Author's Address: University of Setif, The Institute of Optics and Precision Mechanics, 19000 Setif, Algeria, charifi_mohamed@yahoo.fr 103 Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 [6] determined the boundary conditions in the process of binary alloy solidification when the temperature measurements at selected points of the cast are known. In this model of Stefan, the liquidus temperature varies with the concentration of the alloy component. Various direct methods were also used by some authors for determining the temperature at a fixed boundary, knowing the position of a mobile interface solid/liquid [7] or for evaluating the position of one or many mobile interfaces for a given temperature [8] to [10]. Other studies on phase change on materials (PCM) were treated using cylindrical or spherical coordinates [11] to [13], by fixing the mobile front in the presence of convection [14] and increasing energy storage [15]. The position of the solid/liquid interface can be directly determined experimentally using temperature measurements [16], optical monitoring [17], x-ray [18], ultrasound [19] and [20], eddy current [21], thermoelectric [22], and electrical resistance diagnostics. Although these experimental techniques are largely used for identifying the solid/liquid interface in a phase change problem, they remain expensive and difficult to implement. Moravak et al. treated experimentally the solidification conditions of cold worked high alloy tool steel in comparison to construction steels in quasi equilibrium state [23]. Steiner Petrovič and Šturm experimentally study the modification of a non-oriented electrical steel sheet entirely treated with antimony using a laser surface alloy [24]. Hriberšek et al. solved inverse problem for determining surface heat transfer coefficient between liquefied nitrogen and a plate of inconel 718. The design of the numerical simulation was validated experimentally [25]. Indirect identification or control of the solidliquid interface requires a simulation of the problem and its resolution. In the case of material solidification, the phase change problem can be resolved by inverse methods, usually validated by experimental data. From a theoretical perspective, the validation of a simulation (mathematical model used) can be done using exact solutions in particular situations instead of experimental data. In a solidification phase change process, the solid and liquid phases have different thermophysical properties. In the solid phase, the heat transfer is purely conductive. In the liquid phase, two heat transfer situations are possible. If the liquid phase is maintained at the phase change temperature as in our study case, the heat transfer is purely conductive. When the temperature of the liquid phase is greater than the phase change temperature, the heat transfer occurs by diffusion and natural convection. The purpose of this work is to control a pure material solidification by an optimization method formulated as an inverse problem of thermal conduction. We are interested in the control of the solid/liquid interface for a pure material in one-dimensional spherical geometry considering a planar front. The control variable used is the temperature at the fixed boundary of the solid domain. The given problem data are the initial state, the desired front (planar front) evolution, and the phase change temperature. With supplementary information on the heat flux at the front deduced from a heat balance, the inverse problem resolution can be made. The thermal system of our problem is governed by the conduction equation. We introduce for this resolution the least square criteria characterizing the difference between the dynamic behaviour of the system and the corresponding developed mathematical model. We then introduce the adjoint equation for evaluating the criteria gradient. The use of an iterative algorithm based on a conjugated gradient enables to find the optimal solution. We described the problem equation with space-and-time continuous variables and the procedure for finding the criteria gradient with the corresponding adjoint equation. The numerical resolution was undertaken using the conventional finite difference method with a mobile grid related to the considered physical domain. The time variable is then discretized according to the Crank-Nicolson unconditionally stable scheme. The exposed results concern the exact built solution and the noisy exact data. In order to guaranty the well-posed problem and a stable solution for the noisy data case, Tickonov's regularization method was used. 1 INVERSE PROBLEM FORMULATION The equations that govern the problem for the solid and the liquid phases are presented in the following. 1.1 Solid Phase In the solid phase, the temperature T(r, t) re[Sd(t), R] and the heat flux $s(t) penetrating the moving boundary Sd(t) in the time interval [0, f are obtained by the resolution of the following thermal conduction equations: dT(r, t) As d( 2 dT(r, t) PC dt r2 dr I dr T (0, t ) = U (t ), = 0, (1) (2) 104 Charifi, M. - Zegadi, R. Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 T ( S d (t), t ) = T Sd (t = 0) = So, T (r ,0) = T0(r ), dT % (t) = K dT ( s d (t ), t ) dr (3) (4) (5) (6) 1.2 Liquid Phase In the liquid phase, the temperature T(r,t) re[0,Sd(t)] and the heat flux ^(t) penetrating the moving boundary Sd(t) in the time interval [0, tf] are obtained by the resolution of the following thermal conduction equations: dT (r, t) A, d( 2 dT (r, t) Picr dt r2 dr ^ dr T ( S„ (t), t ) = Tf, dT (0, ? ) = 0, dr Sd (t = 0) = S0, T (r ,0) = T0(r ), dT ( (t), t ) = 0, ç,(t) = X, dr (7) (8) (9) (10) (11) (12) where the variables r, t and T in their non-dimensional forms are defined as follows: $ r $ t ■ A r =—, t = - R p-c ■ R2 where Tref = T(0,0) - T(0, tf ). T =- T - T, f (13) ref Remark: thereafter the symbol * will be omitted. During solidification, the heat balance (Stefan equation) at the solid/liquid interface can be expressed as: dS (t ) Vp(t ) (t ) = Vl (t ) = PiL dt (14) In the case in which the liquid phase is maintained at the melting temperature Tf , the heat flux entering the solid is: r dS C) 9p (t)=% (t) = piL~di^~ ■ (15) To solve the direct problem, we need to determine T(r,t;U) and $s(t;U) from the data: {tf,S0, T0(r), Tf,Sd(t)} and from the U(t) control in the interval (0 < t < tf). If ^(t) is the prescribed flux entering the solid at the front, the inverse problem to solved is to find the control U*(t) on [0,tf from the data set in order to obtain $s(t;U*) as close as possible to $p(ty Fixed boundary U(t) Fig. 1. Definition domain solidification process 1.3 Assumptions This simulation concerns a pure material where the following assumptions are usually admitted for this type of material: • The material in each phase is homogeneous and isotropic. • The thermo-physical properties (A,p, c) of the material are independent of temperature, but they are different form phase to phase. • The effect of natural convection in the liquid phase of the material is not taken into account (constant density). • No internal heat generation and all radiation effects are neglected. • The mould wall is considered very thin with no temperature gradient. Its thermal resistance can be neglected. 2 INVERSE RESOLUTION METHOD The introduced least squares criterion J is: J(U) = \'tw1(t)(ys (t;U)-qp (t ))2 dt. (16) The problem resolution consists in determining U*(t) with the conditions: U*eV ,Vt e [0,tf] such that: J(U*) = inf J(U) where V is the set of admissible solutions. A Ribiere-Polak conjugate gradient algorithm is used to calculate the time interval [0, tf] the iterates U(t) according to the following steps: 1. Determination of T (r, t; Un) and $s(t; Un) as solutions of the direct problem in the solid phase. 2. Evaluation of the criterion J( U*). 3. Resolution of the adjoint problem to calculate VJ" (t;U), the gradient of the criterion J with respect to U. 4. Solving the problem of variation to determine 8T(r, t)d. 5. Determination of the iterates Un+1(t) using the following equation: Inverse Method for Controlling Pure Material Solidification in Spherical Geometry 105 Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 Un+1(t) = Un(t) -d"d", n = 1,2,3,... , (17) with the descent direction dn defined as: dn =VJ" (t ;U ) + andn-1, (18) and a" =- JC pJ"] ' d Gateaux's directional derivative DSU(U) of the functional J at point U in the direction SU is defined by: (U) = Dro / (U) = lim J (U + )- J (U ), (19) s^O s SJ(U) = 2j0'/w1(t) (P(t;U)-vp(t)) 8(ps(t;U) dt, (20) DSU J(U) =jVJ(t ;U)8U (t)dt. (21) 0 Using the definition of Gateaux's derivative and computing the variations, we obtain the variation problem where ST and 30 are solutions of the following variation equations: dST(r,t) _ 2 dST(r, t) d2ST(r, t) dt r dr dr2 ST (1, t ) = SU (t ), ST ( Sd (t), t ) = 0, Sd (t = 0) = S0, ST (r ,0) = 0, ÔST(SMJ) 8

k \ 2(Sd -1)A^ (34) To verify the numerical model and the control algorithm reliability, we built the following exact solution: T (r, t) = -r I (0.2)21-(0.2)r J-0.2) (35) 106 (32) Charifi, M. - Zegadi, R. 3.1 Case of No-Noisy Exact Data The optimization algorithm used to determine the optimal control U(t) that ensures a desired evolution front Sd(t) was tested with and without a weighting function. When the weighting function is (wj(t) = 1), the algorithm converges but after a significant number of iterations (n >1000). Instead, the introduction of the weighting function wj(t) = t2 (Fig. 2a) allows the algorithm to converge more quickly after few iterations (n = 54). Consequently, the weighting function has an impact on the computation time and the number of iterations. The improved results over the considered entire thermal process interval depend on the initial guess. The latter must not be chosen arbitrarily. It must belong to the domain of admissible solutions and obey the laws of heat transfer. Figs. 2a and b show that the control U(t) the flux at the fixed boundary are obtained with high precision after a few number of iterations on 95 % of the time horizon. At the remaining time, important errors appear as expected from these global optimization methods. k Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 a) b) Fig. 2. a) Exact control and calculated and b) exact flux and calculated (no-noisy exact data) 3.2 Case of Noisy Exact Data To simulate these errors and generate a data flux ^p(t) close to reality, a white noise b(t) is added to the data ^p(t) with a 5 % maximal amplitude using the following relation: VP (t) =VP (t) exact + b(t) (36) with b(t) = y(t)(0.05)^ref y(t) is a random variable with a uniform probability density over [-1, 1] and $ref is defined as: Vref Vp (tf )exact + Vp (0)e, 2 (37) Fig. 3 shows that the inverse methods are susceptible to errors. The white noise added to the data $p(t) caused significant oscillations on the results. The relative error on the evaluated control U(t) is about 20 %. This is related to the ill-posed nature of the problem. Consequently, these results are not satisfactory. According to [26], in the case of noisy a) b) Fig. 3. a) Exact control and calculated and b) exact flux and calculated (noisy exact data) Inverse Method for Controlling Pure Material Solidification in Spherical Geometry 107 Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 data, the computing stop criterion must satisfy the following condition: | J (Un+1) - J (Un) \2(t ) dU (t) dt dt, (39) where w2(t) = (f-1) is a weighting function and n is a regularization parameter. After derivation, we obtain the following discreet gradient: 2 a) b) Fig. 4. a) Exact control and calculated and b) exact flux and calculated (Regularized problem) O noisy data —*— regularized problem E(t)= 1-U(t)calculated/U(t)exact I V a) 0.2 0.4 0.6 0.6 time (t) [adimensional] b) Fig. 5. a) Relative error on the control, and b) relative error on the flux 108 Charifi, M. - Zegadi, R. Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 VJ(t;U) = -dM--2n [U(t)-w2(t)U"(t)]. (40) dr The choice of parameter n is made heuristically by successive tests while inspecting the solution regularity. The results obtained using Tikhonov's method are satisfactory as shown in Fig. 4. Fig. 5 represents the comparison of the relative error control Fig. 5a and flux Fig. 5b for no-noisy exact data, noisy exact data and those of the regularized problem. The oscillations amplitudes clearly decrease for the regularized problem. The control U (t) and the flux $ (0,t) are stabilized over time. 4 CONCLUSIONS An inverse method to control the evolution of the solidification front of a pure material in a one-dimensional spherical geometry problem was studied. This "ill-posed" problem was solved by using information as a prescribed flux at the front deduced from a heat balance. It was formulated as an optimization problem in which a criterion of least squares is introduced between the model and the object. We then introduced the adjoint equation to accurately calculate the gradient criterion. The obtained equations were discretized using the finite differences classic method in a moving mesh domain. The numerical model was validated by an exact built solution. The resolution was carried out using the global conjugate gradient method. In the case of non-nosy exact data, the algorithm has enabled the rapid control determination with excellent accuracy except for the time steps near the end. The introduction of a selected initial guess has allowed the algorithm to converge more quickly. In the presence of random errors (noisy exact data case) on the prescribed flux, the used regularization method reduced their effects on the obtained control. These obtained results constitute a good starting point for future work in which we consider the pure material solidification resolution in two-dimensional geometry and introducing convection. 5 NOMENCLATURE C specific heat [Jkg-1k-1], d direction of displacement, J criteria, L latent heat [Jkg-1], T temperature [K], P adjoint variable, U control, S position front, V set of admissible solutions, W\, w2 weighting functions, R radial coordinate [m], t time coordinate [s], R dimensionless radial coordinate. Greek symbols a, s positive real numbers, £ transform space coordinate, n regularization parameter, X thermal conductivity [Wm-1 K-1], p material density [kgm-3], d, $ differential, heat flux. Subscripts 0, f initial final 1, k space time nodes index 5, l solid, liquid n iteration number 6 REFERENCES [1] Ozisik, M.N., Orlande, H.R.B. (2000). Inverse Heat Transfer: Fundamentals and Applications. Taylor and Francis, New York. [2] Zabaras, N., Ruan, Y. (1989). A deforming finite element method analysis of inverse Stefan problems. International Journal for Numerical Methods in Engineering, vol.28, no. 2, p. 295-313, DOI:10.1002/nme.1620280205. [3] Zabaras, N., Kang, S. (1993). On the solution of an ill-posed design solidification problem using minimization techniques in finite and infinite-dimensional function spaces. International Journal for Numerical Methods in Engineering, vol. 36, no. 23, p. 3975-3990, DOI:10.1002/nme.1620362304. [4] Samai, M., Jarny, Y., Delaunay, D. (1993). An optimisation method using an adjoint equation to identify solidification front location. Numerical Heat Transfer, Part B: Fundamentals, vol. 23, no. 1, p. 67-89, DOI:10.1080/10407799308914890. [5] Jiang, H., Nguyen, T.H., Prud'homme, M. (2005). Control of the boundary heat flux during the heating process of a solid material. International Communications in Heat and Mass Transfer, vol. 32, no. 6, p. 728-738, DOI:10.1016/j. icheatmasstransfer.2004.10.009. [6] Hetmaniok, E., Stota, D. (2012). Numerical procedure of solving some inverse problem in solidification of the binary alloy. Computer Assisted Methods in Engineering and Science, vol. 19, p. 393-402. [7] Ren, H.-S. (2007). Application of the heat-balance integral to an inverse Stefan problem. International Journal of Thermal Sciences, vol. 46, no. 2, p. 118-127, DOI:10.1016/j. ijthermalsci.2006.04.013. [8] Momose, K., Yamakawa, T.,Kimito, H.(1998). An inverse analysis of two-phase Stefan problems using imaginary heat sources. Heat Transfer - Asian Research, vol. 27, no. 3, p. 179-191, DOI:10.1002/(SICI)1520-6556(1998)27:3<179::AID-HJ1>3.0.CO;2-S. Inverse Method for Controlling Pure Material Solidification in Spherical Geometry 109 Strojniški vestnik - Journal of Mechanical Engineering 63(2017)2, 103-110 [9] Yu, X., Nelson, D.J., Vick, B. (1995). Phase change with multiple fronts in cylindrical systems using the boundary element method. Engineering Analysis with Boundary Elements, vol. 16, no. 2, p. 161-170, DOI:10.1016/0955-7997(95)00052-6. [10] Nowak, I., Nowak, A.J., Wrobel, L.C. (2002). Identification of phase change fronts by Bezier splines and BEM. International Journal of Thermal Sciences, vol. 41, no. 6, p.492-499, D0I:10.1016/S1290-0729(02)01342-X. [11] Regin, A.F., Solanki, S.C., Saini, J.S. (2006). Latent heat thermal energy storage using cylindrical capsule: Numerical and experimental investigations. Renewable Energy, vol.31,no. 13, p.2025-2041, D0I:10.1016/j.renene.2005.10.011. [12] Bilir, L., Ilken, Z. (2005). Total solidification time of a liquid phase change material enclosed in cylindrical/spherical containers. Applied Thermal Engineering, vol. 25, no. 10, p. 1488-1502, D0I:10.1016/j.applthermaleng.2004.10.005. [13] Assis, E., Katsman, L., Ziskind, G., Letan, R. (2007). Numerical and experimental study of melting in a spherical shell. International Journal of Heat and Mass Transfer, vol. 50, no. 9-10, p. 1790-1804, D0I:10.1016/j. ijheatmasstransfer.2006.10.007. [14] Ismail, K.A.R., Da Silva, M.G.E. (2003). Numerical solution of the phase change problem around a horizontal cylinder in the presence of natural convection in the melt region. International Journal of Heat and Mass Transfer, vol. 46, no. 10, p. 1791-1799, D0I:10.1016/S0017-9310(02)00487-8. [15] Lamberg, P. (2004). Approximate analytical model for two-phase solidification problem in a finned phase-change material storage. Applied Energy, vol. 77, no. 2, p.131-152, D0I:10.1016/S0306-2619(03)00106-5. [16] Kittl, J.A., Reitano, R., Aziz, M.J., Brunco, D.P, Thompson, M.O. (1993). Time-resolved temperature measurements during rapid solidification of Si-As alloys induced by pulsed-laser melting. Journal of Applied Physics, vol. 73, no. 8, p. 37253733, D0I:10.1063/1.352903. [17] Queheillallt, D.T., Lu, Y., Wadeley, H.N.G. (1997). Laser ultrasonic studies of solid-liquid interface. Journal of the Acoustical Society of America, vol. 101, no. 2, p. 843-853, D0I:10.1121/1.418042. [18] Kakimoto, K., Eguchi, M., Watanabe, H. Hibiya, T. (1988). In-situ observation of solid-liquid interface shape by X-ray radiography during silicon single crystal growth. Journal of Crystal Growth, vol. 91, no. 4, p. 509-514, D0I:10.1016/0022-0248(88)90118-2. [19] Parker, R.L., Manning, J.R., Peterson, N.C. (1985). Application of pulse-echo ultrasonics to locate the solid/liquid interface during solidification and melting of steel and other metals. Journal of AppliedPhysics, vol. 58, no. 11, p. 4150-4164, D0I:10.1063/1.335547. [20] Burhan, D., Ihara, I., Seda, Y. (2005). In situ observations of solidification and melting of aluminum alloy using ultrasonic waveguide sensor. Materials Transactions, vol. 46, no. 9, p. 2107-2113, D0I:10.2320/matertrans.46.2107. [21] Wadley, H.N.G., Dharmasena, K.P. (1997). Method for liquidsolid interface shape and location discrimination during eddy current sensing of Bridgman growth. Journal of Crystal Growth, vol. 172, no. 3-4, p. 313-322, D0I:10.1016/S0022-0248(96)00496-4. [22] Verhoeven, J.D., Gibson, E.D. (1969). A thermoelectric probe for measuring solid-liquid interface velocities and temperatures. Journal of Crystal Growth, vol. 5, no. 4, p. 235238, D0I:10.1016/0022-0248(69)90050-5. [23] Moravcík, R., ¿tefániková, M., Čička, R., Čaplovič, L., Kocúrová, K., Šturm, R. (2012).Phase transformations in high alloy cold work tool steel. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 12, p. 709-715, D0I:10.5545/sv-jme.2012.531. [24] Steiner Petrovič, D., Šturm, R. (2014). Fine-structured morphology of a silicon steel sheet after laser surface alloying of sb powder. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 1, p. 5-11, D0I:10.5545/sv-jme.2013.1347. [25] Hriberšek, M., Šajn, V., Pušavec, F., Rech. J., Kopač, J. (2016). The procedure of solving the inverse problem for determining surface heat transfer coefficient between liquefied nitrogen and Inconel 718 workpiece in cryogenic machining. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 6, p. 331-339, D0I:10.5545/sv-jme.2016.3572. [26] Alifanov, O.M., Mikhailov, V.V. (1978). Solution of the non-linear inverse thermal conductivity by the iteration method, Journal of Engineering Physics, (Translated from Inzhenerno-Fizicheskii Zhurnal, vol. 35, no. 6, p. 1123-1129), D0I:10.1007/BF01104861. [27] Tikhonov, A., Arsenine, V.Y. (1977). Solution of III-Posed Problems, Wiley, New York. 110 Charifi, M. - Zegadi, R.