A method for transient thermal load estimation and its application to identification of aerodynamic heating on atmospheric reentry capsule

Aerospace Science and Technology 51 (2016) 26–33

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

A method for transient thermal load estimation and its application to identification of aerodynamic heating on atmospheric reentry capsule Piotr Duda Cracow University of Technology, Institute of Thermal Power Engineering, Al. Jana Pawła II 37, 31-864 Kraków, Poland

a r t i c l e

i n f o

Article history: Received 9 October 2015 Received in revised form 19 January 2016 Accepted 22 January 2016 Available online 25 January 2016 Keywords: Inverse method Transient heat flow Identification Reentry capsule

a b s t r a c t The purpose of this work is to formulate a simple method which can be used for the transient thermal load identification. The proposed algorithm can reconstruct the transient temperature distribution in a whole component based on measured temperatures in selected points on the component surface. The presented method may use some commercial heat transfer modelling software. In this work, ANSYS Multiphysics is applied. Two numerical examples of identification of the transient thermal load taken from literature are presented. The determination of the heating process of a simple two-dimensional plate and the recognition of the aerodynamic heating on the atmospheric reentry capsule will demonstrate high accuracy and stability of the proposed algorithm. Future time steps and smoothing filters are used to stabilize the solution of an ill-posed problem. The presented method is easy to use for engineers working in the aerospace or power industry. © 2016 Elsevier Masson SAS. All rights reserved.

1. Introduction The thermal analysis plays a significant role in the design of components subjected to high temperatures and intended for many engineering fields such as the aerospace or power industries. One of the major difficulties in carrying out the thermal analysis is the definition of all boundary conditions. The boundary condition formulation is not easy due to many factors, e.g. the irregular geometry of the component edge or the trouble in interpreting the physical phenomena occurring on the surface in time and space. An interesting method of solving this difficulty is to include additional temperature measurements instead of unknown boundary conditions and then formulate a transient inverse heat conduction problem (IHCP). For elements with simple and regular shapes, assuming that the material thermal and physical properties are constant, exact methods of solving the 1-D inverse heat conduction problem may be used [1–3]. Solving multidimensional inverse problems related to bodies with temperature-dependent physical properties requires numerical methods. The methods of solving steady- and transientstate two-dimensional inverse problems of heat conduction in solid bodies with simple shapes are presented in [4] and [5–7], respectively. The algorithms used for inverse problems in three dimensions are presented in [8,9] and a general method for solving multidimensional inverse heat conduction is shown in [10]. There are very few works presenting the solution of multidimensional

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.ast.2016.01.015 1270-9638/© 2016 Elsevier Masson SAS. All rights reserved.

inverse heat conduction in solid bodies with complex shapes. An attempt is made in [11] to use the method [12] to identify the unknown heat flux on the surface of a reentry capsule with a complex shape, but the number of necessary temperature measurements is very high and difficult to carry out in practice. Inverse methods may be also used for optimization of the power unit start-up and shut-down operations, which may lead to the heat loss reduction during these operations and extend the power unit life span [13,14]. The IHCP solution often includes a complex mathematical model [1–14] and therefore it is difficult to perform for many engineers. Commercial heat transfer modelling software is hard to apply because this requires that all boundary conditions should be defined. The purpose of this work is to propose a simple inverse method which can be used for the transient thermal load estimation on solid bodies with complex shapes. The biggest advantage of the proposed method in comparison with the existing ones is its simplicity. This original new formulation does not include a sophisticated mathematical model and may use some commercial heat transfer modelling software. In this work, ANSYS Multiphysics software [16] that utilizes the finite element method (FEM) is applied. Another advantage is the generality of the proposed formulation. This algorithm will be applied to identify aerodynamic heating on the atmospheric reentry capsule, but it may also be used for many other transient phenomena if they have a direct solution only. The method may be used to solve multidimensional inverse problems in bodies with simple and complex geometries. It makes it possible to take account

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

of different nonlinearities, such as temperature-dependent material properties or nonlinear boundary conditions, and this paper presents radiation-related nonlinearity as well. In this work, future time steps [15] and smoothing filters [14] are used to stabilize the ill-posed problem solution. Two numerical examples of identification of the transient thermal load taken from [11] will be presented. The determination of the heating process of a simple two-dimensional plate and the recognition of the aerodynamic heating on the atmospheric reentry capsule will demonstrate high accuracy and stability of the proposed algorithm. The results obtained by the proposed method and the method presented in [11] will be compared.

27

Fig. 1. A direct transient heat conduction problem.

2. Formulation of the method The equation governing the transient heat conduction in solids is expressed as follows:

c ( T )ρ ( T )

∂T = −∇ · q, ∂t

(1)

where q is the heat flux vector defined by Fourier’s law

q = −D∇ T .

(2)

D is the conductivity matrix

 D=

kx ( T ) 0 0 0 k y (T ) 0 0 0 kz (T )

 (3)

If the material is isotropic, then kx ( T ) = k y ( T ) = k z ( T ) = k( T ). All the material properties: c – specific heat, k – thermal conductivity, ρ – density can be considered as known temperature or temperature-independent functions. Transient heat conduction problems are initial-boundary problems and appropriate initial and boundary conditions have to be defined. The initial condition is the temperature of a body at its first moment t 0 = 0 s.

T (r , t )|t0 =0 = T 0 (r )

Fig. 2. An inverse problem with known and unknown boundary conditions, and additional temperature measurement points.

(4)

where r is the position vector. Typically, three boundary conditions – of the 1st, 2nd and 3rd kind – may be assigned to the body boundary

T | T = T b

(5)

( D ∇ T · n)|q = q B

(6)

( D ∇ T · n)|h = h( T m − T |h )

(7) If boundary conditions (5)–(8) are unknown on parts of the domain boundary, the problem becomes ill-posed and additional interior temperature measurements are needed in the analysis (Fig. 2)

where n – unit outward normal vector to boundary  , T b – temperature set on the body boundary  T , q B – heat flux set on the body boundary q , h – heat transfer coefficient set on the body boundary h , T m – temperature of the medium.

f i (t ) = T (r i )

If a body surface is cooled or heated by radiation, then the boundary condition is nonlinear

  ( D ∇ T · n)|r = σ F T r4 − ( T |r )4

Fig. 3. Staircase or piecewise linear function of an unknown heat flux on the body surface.

(8)

where σ is a Stefan–Boltzmann constant, T r is the surrounding temperature and F is a shape coefficient which takes account of the surface emissivity and the mutual configuration function. A phenomenon described by equations (1)–(8), in which boundary conditions are specified over the entire boundary, is referred to as a direct transient heat conduction problem (Fig. 1).

i = 1, . . . , N T

(9)

The first step of the proposed method is the discretization of any unknown boundary condition type into N u intervals and its approximation by the 2nd kind boundary condition (6). It may be assumed as a staircase or a piecewise linear function on the body surface (Fig. 3). The aim is to choose such u (u 1 (t ), u 2 (t ), . . . , u N u (t )) in time that the computed temperatures should agree within certain limits with the temperature values measured experimentally. It may be expressed as

T i (u , r i , t ) − f i (t ) ∼ = 0,

i = 1, . . . , N T

(10)

28

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

The Levenberg–Marquardt method (LM) is a combination of the Gauss–Newton method (μ(k) → 0) and the steepest-descent method (μ(k) → ∞). The iterative procedure is continued until (k)

changes in u L are less than  . At every k-th iteration step, the time-invariant vector u L is calculated, but temperature distribu(k) tion T (u L , ri , t ) is time-dependent. The initial-boundary value problem given by Eqs. (1)–(8) is solved at each iteration step by the finite element method (FEM) calculating qB from u L . The computational procedure for the solution of the inverse problem using the LM method is summarized as follows:

Fig. 4. Diagram illustrating calculation of u L using future time steps.

The number of measurement points N T should not be smaller than the number of unknown boundary conditions N u . Ill-posed problems are hard to solve. They are very sensitive to random measuring errors in input data. Many methods may be used to stabilize the solution. In this paper, future time steps and smoothing filters are used. In order to stabilize the solution, vector u L is kept constant in N F future time steps (Fig. 4)

u L = u L +1 = · · · = u L + N F −1

(11)

where u L = u (u 1 (t L ), u 2 (t L ), . . . , u N u (t L )). The application of this well-known algorithm (future time steps) to stabilize the solution makes it possible to formulate the proposed inverse method. The least-squares method is used to determine parameters u L . The sum of squares

S=

N T L + N F −1   i =1

2

f (t j )i − T (u L , r i , t j )

=

j=L

NT 

F i (u L )

(12)

i =1

may be minimized by a general unconstrained method. However, the properties of (12) make it worthwhile to use methods designed specifically for the nonlinear least-squares problem. In this paper, the Levenberg–Marquardt method [17] is used to determine parameters u L . Application of the Levenberg– Marquardt method for the identification of thermal boundary conditions in water-wall tubes of boiler furnaces in steady state conditions can be found in [21]. It is assumed that the searched-for vector changes in the range

ulL ≤ u L ≤ u uL

(13)

where superscripts l and u denote the lower and upper bounds of the searched-for parameters, respectively. The Levenberg–Marquardt method performs the kth iteration as (k+1)

uL

(k)

= u L + δ (k) ,

(14)

where

δ

(k)

=



J

 (k) T

J

(k)

(k)

+μ

In

−1 

J

 (k) T



(k) 

F uL

k = 0, 1 , . . . ,

, (15)

where μ is the multiplier, In is the identity matrix and J is the Jacobian matrix, which has the following form

⎡ ⎢ ⎢ ⎢ (k) J =⎢ ⎢ ⎣

∂ F1 ∂ u 1 (t L )

... ... ... ...

∂ F NT ∂ u 1 (t L )

... ... ... ... ... ...

∂ F1 ∂ u N u (t L )

... ... ... ...

∂ F NT ∂ u N u (t L )

⎤   ⎥ ⎥ ⎥ ⎥ ⎥ ⎦  

1. Solve the initial-boundary problem defined by Eqs. (1)–(8) in the time interval t ∈ 0, t N F assuming the starting value of the searched-for vector u L ( L = 0) from the range defined by (13). 2. Calculate the objective function vector F (u L ( L = 0)) defined in (12). 3. Use the LM method to find vector u L ( L = 0) which minimizes the objective function vector F (u L ( L = 0)). 4. Save the calculated vector u L ( L = 0) and the temperature distribution at time t = t N F . 5. Solve the initial-boundary problem defined by Eqs. (1)–(8) in the time interval t ∈ t N F , 2t N F assuming the starting value of vector u L ( L = N F ) = u L ( L = 0) and set the temperature distribution calculated at time t = t N F for the initial condition (4). 6. Use the LM method to find vector u L ( L = NF ) which minimizes the objective function vector F (u L ( L = N F )). 7. Save vector u L ( L = NF ) and the temperature distribution at time t = 2t N F . 8. Continue the procedure to the end of the transient process (t end − t N F ). The method presented for the solution of the inverse heat conduction problem may be generalized to solve different inverse transient heat transfer problems. The applicability of the presented algorithm will be demonstrated in the solution of inverse heat conduction and simultaneous inverse heat transfer by conduction and radiation for the reentry capsule. 3. A square plate heating identification The square plate ABCD presented in Fig. 5 with the dimension L = A B = 0.5 m is subject to the following boundary conditions [11]

q(x, y , t ) =

⎧ ⎪ ⎨ q (t ),

on AA’



− q (t ), on CC’ ⎪ ⎩ on other boundary 0,

and

⎧ 0, ⎪ ⎪ ⎨ Q (t − t )/(t − t ),  1 0 1 0 q (t ) = ⎪ Q , ⎪ ⎩ 1 Q 2,

(k)

u L =u L

0 ≤ t ≤ t0 t0 ≤ t ≤ t1 t1 ≤ t ≤ t2

(18)

t2 ≤ t ≤ t3

where the coordinates of A, A , C , C ’ are A (0, 0), A ( L /2, 0), C ( L /2, L ), C ( L , L ), respectively, and L = 0.5 m, t 0 = 50 s, t 1 = 100 s, t 2 = 300 s, t 3 = 450 s, Q 1 = 50 kW/m2 , Q 2 = 100 kW/m2 . 

(16)

(17)

The heat flux history q (t ) is presented in Fig. 7 using asterisks. The component is made of titanium, which has the following material properties: thermal conductivity k is 21.9 W/mK, specific heat c is 520 J/kg K and density ρ is 4507 kg/m3 . Initially, the plate has a uniform temperature T 0 = 0 ◦ C.

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

Fig. 5. Discretization of a square plate used in the direct numerical (FEM) and inverse solution.

29

transients and two given measured temperature transients. Vector u L includes two components u 1 (t L ), u 2 (t L ). The proposed method is tested using the “smoothed noisy measured data” and the time step t = 1 s. For the lower and upper bounds of all searched-for parameters, the values 0.1 and 150 kW/m2 are chosen in vector u L . The starting vector u L ( L = 0) = 2 kW/m2 , 2 kW/m2 is assumed. Additionally, the number of future time steps N F = 10 is selected. The computational procedure presented in chapter 2 was programmed in the Fortran 90 language. Microsoft IMSL Libraries were used to perform the Levenberg–Marquardt algorithm [20]. The Jacobian matrix defined by (16) was calculated in this subroutine automatically, based on the objective function vector F defined in (12). ANSYS Multiphysics software was called as a subroutine from the main program. The number of iterations at each time step varies between 1 and 48. The calculated heat flux histories as the staircase function can be seen in Fig. 7. If the midpoints of the constant heat flux intervals are connected using straight lines, the identified boundary condition is close to the exact one. This method of presenting the calculated heat flux history will be used in this paper. For a quantitative accuracy estimation of the proposed inverse method, the following error measure is defined

E=

N time N u 

  inv u (t j ) − u exact (t j ) i

j

i

(19)

i =1 t

Fig. 6. “Noisy measured data” and “smoothed noisy measured data” calculated using the FEM.

The analysed domain is divided into 64 linear plane finite elements and t = 1 s is selected as the time step. The transient temperature distribution is calculated in the whole plate using the ANSYS Multiphysics software [16]. Temperature transients f 1 (t ) and f 2 (t ) are assumed as “exact measured data” for the inverse solution. In order to bring the numerical test closer to real conditions, “noisy measured data” are obtained by adding normal random errors ±2 ◦ C of the zero mean to the “exact measured data” (σf = 2/3), and they are presented in Fig. 6. An ill-posed problem is difficult to solve, as calculated temperatures are very sensitive to measurement errors. The errors may create temperature oscillations, which may lead to an unstable solution. Gram polynomials are used to rectify temperature data [14]. The value of the Gram smoothing polynomial is calculated at the centre of the approximation interval, while for the determination of polynomial coefficients – future and past data points are used. In order to reduce the oscillations, the “noisy measured data” are smoothed out using the above-mentioned filters based on the local polynomial approximation (3rd order and eleven subsequent measurement points). The inverse problem is then formulated without specifying boundary conditions on edges AA’ and C’C, while two measurements f 1 (t ) and f 2 (t ) are added (see Fig. 5). It is assumed that the unknown boundary condition on edge AA’ consists of two unknown heat fluxes u 1 (t ), u 2 (t ). The unknown boundary condition on edge C’C consists of −u 2 (t ), −u 1 (t ). There are two unknown

where N time = tend , u inv (t j ) – heat flux calculated using the proi NF posed inverse method, u exact (t j ) – heat flux obtained from the i direct solution. Error E of the identified heat flux transients for the “smoothed noisy measured data” equals 3.52 kW/m2 . If the absolute error E is related to the maximum heat flux of 100 kW/m2 during the heating process, the relative error may be calculated – it equals 3.52%. Smoothing filters and future steps are very effective in eliminating random measurement errors. The use of Gram polynomials to smooth the “noisy measured data” does not allow a stable calculation by means of the proposed method. It is necessary to use also a suitable number of future time steps N F . The first stable solution is obtained for N F = 9 and the time step t = 1 s. An excessive increase in the N F number involves a rise in the solution errors. The lowest error is observed for N F = 10 and the time step t = 1 s. In general, the optimum value of N F can be found by performing a series of calculations for a different number of future time steps N F and selecting the smallest number that ensures stabilization of the solution. Fig. 8 presents a comparison of unknown boundary conditions estimated by the proposed method with a different number of future time steps N F . During each iteration the initial-boundary problem defined by Eqs. (1)–(8) is solved by the FEM using the ANSYS Multiphysics software. This parametric formulation of the method of solving the transient inverse heat conduction problem simplifies its use. The proposed algorithm may be applied to all transient nonlinear problems if they have a direct solution. 4. Reentry capsule The Japan Aerospace Exploration Agency (JAXA) conducted a conceptual study of a reentry testbed “PARTT” (Piggyback Atmospheric Reentry Technology Testbed). The weight and the diameter of the PARTT reentry module were designed as 40 kg and 0.7 m, respectively. The reentry module materials and geometry are presented in Fig. 9 [11]. The reference heat flux transient of a 0.7 m diameter sphere and the heat flux distribution on the shell are taken from [11] and

30

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

Fig. 7. A comparison between exact heat fluxes and those calculated using the inverse method based on “smoothed noisy measured data”; (a) – heat flux u 1 , (b) – heat flux u 2 .

Fig. 8. Unknown boundary condition estimated by an inverse method with a different number of future time steps N F using the “smoothed noisy measured data”; (a) – heat flux u 1 , (b) – heat flux u 2 .

Fig. 9. Reentry module structure.

presented in Fig. 10a and Fig. 10b, respectively. For further analysis, the spatial distribution of the heat flux is also approximated using a piecewise linear function. The heat transfer between the hot structure and the cold air around the capsule is ignored [11].

The transient temperature distribution in the reentry module structure for the inverse analysis was prepared by a direct FE analysis with the heat flux presented in Fig. 10. The FE model is composed of 473 quadratic quadrilateral elements with eight nodes to improve the solution accuracy [18]. Fig. 11 shows a detailed FE model. The material properties are assumed as temperatureindependent as in [11]. Initially, the reentry capsule has a uniform temperature T 0 = 20.5 ◦ C, and t = 1 s is chosen for the time step in the transient analysis. The problem is nonlinear due to the radiative heat transfer on the “known boundaries” of the model. Temperature histories in selected points and the temperature distribution after 210 s of the heating process can be seen in Fig. 12 and Fig. 13, respectively. The obtained results are similar to the values presented in [11,19]. The proposed inverse method is used for the transient thermal load estimation on the outer surface of the reentry module shell. The shell is presented between points A and C in Fig. 11 and it can also be seen in Fig. 9. The solution is possible thanks to “measured” temperatures in points f 1 , f 2 , f 3 , f 4 , whose location is marked in Fig. 11. The “noisy measured data” are obtained by

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

Fig. 10. Heat flux on the reentry shell; (a) reference heat flux history, (b) distribution on the shell.

Fig. 11. Axisymmetric FE model and location of measurement points f 1 , f 2 , f 3 , f 4 .

adding normal random errors ±1 ◦ C of the zero mean to the “exact measured data” (σf = 1/3) and they are presented in Fig. 14. In order to reduce oscillations, the “noisy measured data” are smoothed out using digital filters. To show the advantages of the presented method, the algorithm will be applied in two ways.

Fig. 12. Temperature of the reentry module structure (direct FE analysis).

4.1. First case It is assumed that the heat flux spatial distribution is known (Fig. 10b), but the reference heat flux transient shown in Fig. 10a is not. The problem can be solved by adding two measured temperature transients f 1 and f 2 . Now, there is one unknown transient and two given measured temperature transients. Vector u L consists of one component u 1 (t L ) specifying the reference heat flux history. For the lower and upper bounds of the searched-for parameter, in vector u L , the values of 0 and 1000 kW/m2 are selected. The starting vector is assumed as u L ( L = 0) = 50 kW/m2 . Additionally, the number of future time steps N F = 20 is chosen. The number of iterations at each time step varies between 1 and 21. The calculated heat flux history can be seen in Fig. 15. The identified reference heat flux transient is close to the exact one.

Fig. 13. Temperature distribution at 210th s (T in K).

31

32

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

Fig. 14. “Noisy measured data” and “smoothed noisy measured data” for the reentry module calculated using the FEM.

Fig. 16. Heat flux on the reentry shell; a comparison between the exact transient and the history calculated using the inverse method based on the “smoothed noisy measured data”.

Fig. 15. Reference heat flux, a comparison between the exact transient and the history calculated using the inverse method based on the “smoothed noisy measured data”.

Fig. 17. Temperature of the reentry module structure; a comparison between the exact temperature transients and the temperature transients calculated using the inverse method.

4.2. Second case

The calculated heat flux transients can be seen in Fig. 16. The identified heat flux transients are close to the histories obtained from the direct solution. Some instability appears only for u 1 after about 300 s of the heating operation. Based on the identified boundary condition, the temperature distribution in the reentry module structure is calculated. The temperature transients obtained from the proposed inverse method are presented in Fig. 17. The reconstructed transient temperature distribution resulting from the inverse method is verified by a comparison with the temperature transients obtained from the direct method. Larger differences appear only in point A at the time of about 300 s. They are caused by the solution instability, which can be seen in Fig. 16.

It is assumed that both the heat flux spatial distribution and the heat flux transient at the stagnation point are not known. The unknown 2nd-kind boundary condition on the outer surface of the reentry module shell is approximated using a piecewise linear function, as presented in Fig. 10b. The problem is solved by adding four measured temperature transients f 1 , f 2 , f 3 and f 4 . Now, there are three unknown transients and four given measured temperature transients. Vector u L consists of three components u 1 (t L ), u 2 (t L ), u 3 (t L ) specifying the heat flux transients at r = 0 mm, r = 325.5 mm and r = 350 mm. For the lower and upper bounds of all searched-for parameters, in vector u L , the values of 0 and 1000 kW/m2 are selected. The starting vector is assumed as u L ( L = 0) = 30 kW/m2 , 40 kW/m2 , 25 kW/m2 . Additionally, the number of future time steps N F = 20 is chosen. The number of iterations at each time step varies between 1 and 45.

5. Conclusion The inverse method presented herein allows identification of the transient temperature distribution and the local heat flux on unknown boundary edges. The unknown boundary condition may

P. Duda / Aerospace Science and Technology 51 (2016) 26–33

be assumed as a staircase function or a piecewise linear function. Introduction of a linear function makes it possible to reduce the number of sought parameters. The method may be applied in practice in an off-line mode. The method can calculate the temperature distribution based on temperature measurements which may be located on easily accessible surfaces. The number of the temperature measurement points is smaller than in the method presented in [11]. It is partly due to the small number of the degree of freedom for the spatial distribution (3 for the present case and 55 for [11]). The proposed algorithm allows the use of commercial calculation programs. In this work, the ANSYS Multiphysics software that utilizes the finite element method (FEM) is applied. Both linear and nonlinear finite elements are used. Due to the implementation of the finite element method, the algorithm may be implemented in respect of complex-shape bodies with temperature-dependent thermophysical properties. The presented inverse procedure is a lot simpler than the method proposed in [11]. The method is stable in terms of data obtained from measurements which are often disturbed by measurement errors. The algorithm applicability is demonstrated by comparing it with the results obtained from the direct and the inverse solution. The presented examples prove good accuracy and stability of the proposed method. Based on the identified transient temperature distribution in the equipment used in the aerospace industry and in power engineering, the transient operation of selected components may be changed to reduce the maximum temperature differences and thermal stresses. As a result, the proposed method may extend the life span of many components. Conflict of interest statement

[3]

[4] [5] [6] [7]

[8]

[9]

[10]

[11]

[12] [13]

[14]

[15] [16] [17] [18]

There is no conflict of interest. [19]

References [20] [1] O.R. Burggraf, An exact solution of the inverse problem in heat conduction theory and applications, ASME J. Heat Transf. 86 (August 1964) 373–382. [2] D. Langford, New analytic solution of the one-dimensional heat equation for temperature and heat flow rate both prescribed at the same fixed boundary

[21]

33

(with applications to the phase change problem), Q. Appl. Math. 24 (4) (1967) 315–322. J. Stefan, Über die Theorie der Eisbildung insbesondere über die Eis-bildung im Polarmeere, Sitzungsber. Kaiserlich. Akad. Wiss. Wien, Math. Naturwiss. Kl 98 (2a) (1889) 965–973. E. Hensel, R. Hills, Steady-state two-dimensional inverse heat conduction, Numer. Heat Transf., Part B 15 (1989) 227–240. H.R. Busby, D.M. Trujillo, Numerical solution to a two-dimensional inverse heat conduction problem, Int. J. Numer. Methods Eng. 21 (1985) 349–359. D.M. Trujillo, Application of dynamic programming to the general inverse problem, Int. J. Numer. Methods Eng. 12 (1978) 613–624. R. Dou, Z. Wen, G. Zhou, 2D axisymmetric transient inverse heat conduction analysis of air jet impinging on stainless steel plate with finite thickness, Appl. Therm. Eng. 93 (2016) 468–475. O.M. Alifanov, A.V. Nenarokomov, Three-dimensional boundary inverse heat conduction problem for regular coordinate systems, Inverse Probl. Eng. 7 (4) (1999) 335–362. S. Lu, Y. Heng, A. Mhamdi, A robust and fast algorithm for three-dimensional transient inverse heat conduction problems, Int. J. Heat Mass Transf. 55 (2012) 7865–7872. Y. Jarny, M.N. Ozisik, J.P. Bardon, A general optimization method using adjoint equation for solving multidimensional inverse heat conduction, Int. J. Heat Mass Transf. 34 (11) (1991) 2911–2919. T. Nakamura, Y. Kamimura, H. Igawa, Y. Morino, Inverse analysis for transient thermal load identification and application to aerodynamic heating on atmospheric reentry capsule, Aerosp. Sci. Technol. 38 (2014) 48–55. P. Duda, Solution of multidimensional inverse heat conduction problem, Heat Mass Transf. 40 (2003) 115–122. A. Rusin, G. Nowak, M. Lipka, Practical algorithms for online thermal stress calculations and heating process control, J. Therm. Stresses 37 (11) (2 November 2014) 1286–1301. P. Duda, Numerical and experimental verification of two methods for solving an inverse heat conduction problem, Int. J. Heat Mass Transf. 84 (2015) 1101–1112. J.V. Beck, Nonlinear estimation applied to the nonlinear inverse heat conduction problem, Int. J. Heat Mass Transf. 13 (1970) 703–716. ANSYS User’s Manual, Revision 14.0 A. G.A.F. Seber, C.J. Wild, Nonlinear Regression, Wiley, New York, 1989. P. Duda, Finite Element Method formulation in polar coordinates for transient heat conduction problems, J. Therm. Sci. 25 (2) (2016), accepted for publication. T. Nakamura, K. Fujii, Probabilistic transient thermal analysis of an atmospheric reentry vehicle structure, Aerosp. Sci. Technol. 10 (2006) 346–354. User’s Manual: IMSL MATH/LIBRARY: FORTRAN Subroutines for Mathematical Applications, IMSL, Houston, Tex., 1989. Print. P. Duda, J. Taler, A new method for identification of thermal boundary conditions in water-wall tubes of boiler furnaces, Int. J. Heat Mass Transf. 52 (2009) 1517–1524.