Determination of transient temperature fields in thick-walled elements using the inverse method

The purpose of this paper is to propose a method of determining the transient temperature of the inner surface of thick-walled elements. The method can be used to determine thermal stresses in pressure elements. An inverse marching method is proposed to determine the transient temperature of the thickwalled element inner surface with high accuracy. The inverse method was validated computationally. The comparison between the temperatures obtained from the direct heat conduction problem solution and the results obtained by means of the proposed inverse method is very satisfactory. The advantage of the method is the possibility of determining the heat transfer coefficient at a point on the exposed surface based on the local temperature distribution measured on the insulated outer surface. The heat transfer coefficient determined experimentally can be used to calculate thermal stresses in elements with a complex shape. The proposed method can be used in online computer systems to monitor temperature and thermal stresses in thick-walled pressure components because the computing time is very short.


Introduction
Unsteady-state thermal stresses arising in pressure elements cannot be calculated correctly without knowing the temperature of the fluid and the heat transfer coefficient on the inner surface of these elements.In order to determine the heat transfer coefficient, very precise measurements of the transient temperature are required due to the small difference between the temperature of the fluid under high pressure and the temperature of the element inner surface.The heat transfer coefficient has a significant effect on the optimum rate of changes in the fluid temperature, which is determined from the condition of not exceeding the thermal stress value on the inner surface of the pressure element [1].
The temperature measurement of the pressure element surface which is in contact with the fluid is difficult, especially if the fluid temperature, pressure and velocity parameters are high.It is not always possible to mount a temperature sensor on the element surface, and temperature measurements can be disturbed by significant errors due to contact resistance.
The most common method for determining timedependent changes in the element surface temperature described in reference literature is based on solving the direct heat conduction problem analytically or numerically.This method is based on the heat transfer coefficient, determined from the correlation for the Nusselt number.The analysis usually concerns one-dimensional cases.Two-dimensional or threedimensional issues are discussed in [2][3][4].Using the direct method to determine the inner surface transient temperature can be problematic.The heat transfer coefficients on the inner surfaces of elements are determined from correlations known from reference literature.These are usually correlations found for hydrodynamically and thermally developed flows under steady-state conditions.
In recent years, extensive research has been conducted on the development of methods based on solving inverse heat conduction problems [5][6][7][8][9][10][11][12][13].These methods enable determination of the inner surface temperature on the basis of two known boundary conditions on the outer surface.In the case of elements with regular shapes, the two-or three-dimensional transient temperature field can be determined by solving the inverse heat conduction problem using marching methods [7,8,11,12].In order to calculate the inner surface transient-state temperature, the inverse region is divided into control volumes.The method marches in space towards the element inner surface by using the energy balance equations for the finite volumes placed on the boundary to find the temperature values in adjacent nodes.
A modified space marching method for transientstate nonlinear inverse heat conduction problems is presented in [9].It is based on the classical space marching method and reduces numerical oscillation by using future temperatures in the least square manner.However, for elements with complex shapes and with holes, identification of thermal stresses and temperature fields by means of the marching inverse method is difficult.In this article, a new approach to temperature and thermal stress determination is proposed.Using the heat transfer coefficient determined on the element inner surface and the fluid temperature [12], commercial programs based on the finite element method can be applied to calculate thermal stresses.In this way, thermal stresses on the edges of holes or other places of stress concentration can be found.

Identification of the temperature distribution in the flat thick-walled plate
An inverse space marching method is used to determine the temperature distribution inside the plate (cf. Figure 1).The transient heat conduction equation for a flat plate is where T denotes temperature, c, ρ, k − specific heat, density and thermal conductivity of the plate material, respectively, and v q & − energy generation rate per unit volume.
the temperature distribution inside the plate including the inner surface (z = s) is found.In the developed calculation method it is assumed that the outer surface, where the temperature is measured, is thermally insulated.Therefore, the boundary condition for z = 0 is as follows: There are two known boundary conditions on the surface z = 0, i.e. conditions ( 2) and (3), while the boundary condition on the opposite surface of the plate z = s is unknown.
The side surfaces of the plate are thermally insulated, so the boundary conditions of the second type are as follows: 1 0 0, 0 The solution of the inverse problem ( 1)-( 5) is obtained using the finite volume method (FVM).
The plate division into finite volumes is shown in Figure 2. The energy balance equations can be written using the control volume method for all nodes located in the centres of the control volumes.

Fig. 2. Plate division into control volumes
The energy balance equation for node (i, j, k) inside the analysed finite volume is: , In Eq. ( 6) the symbols Δx, Δy, Δz denote the control volume dimensions, c(Ti,j,k), ρ(Ti,j,k), k(Ti,j,k) − specific heat, density and thermal conductivity in temperature Ti,j,k, respectively, and , , ( ) x y z q & − energy generation rate per unit volume at points (xi, yj, zk).
Figure 3 shows examples of nodes together with the surrounding cells, for which the heat balance equations are written.Figure 3(a) illustrates all finite volumes taken into account when the heat balance equation is written for a node situated on an insulated surface (z = 0).The thickness of control volumes adherent directly to the insulated surface is Δz/2.Figure 3  ( ) In order to solve the inverse problem, the analysed thickwalled plate is divided into finite volumes (Figures 2-5).First, to solve the inverse problem, the heat balance equations are written for nodes located on the insulated surface z = 0 (Figures 3(a  Using the calculation procedure described above, only the temperature at point (4,4,4) is found.In a similar way, the calculation procedure (calculation template) can be used to determine the temperature of additional nodes lying in plane z = 3Δz = s.In order to determine the heat flux density at point (4,4,4), it is necessary to measure the temperature in 16 additional nodes lying on the insulated surface (z = 0).In this case, it is possible to determined temperature values in five nodes on the plate surface z = 3Δz = s: (4,4,4), (3,4,4), (5,4,4), (4,3,4), and (4,5,4).The flux density of the heat transferred to the fluid from the plate surface in node (4,4,4) can then be found from the heat balance equation for node (4,4,4).The heat transfer coefficient h on surface z = 3Δz = s is determined from the following formula: ( ) , , ( , , ) where Tf denotes the fluid temperature.

Computational validation of the inverse method
The calculations are carried out for a flat plate.It is assumed that the nodes on both surfaces (z = 0 and z = s = 0.045 m) and inside the plate are arranged as in Figures 3 and 4 In the calculations, it is assumed that the density of the heat flux transferred through the plate is q & = 20 kW/m 2 , and the initial temperature T0 is 20ºC.At first, to validate the developed method, "measured data" are generated using the direct method, assuming uniform heating of the plate inner surface (z = s).

Description of the direct method
The "measurement data" are obtained using an analytical method.Because the plate inner surface heating is uniform, the heat conduction equation for the flat wall with independent physical properties of the plate material is where 0 ≤ z ≤ s.Symbol κ = k/(cρ) in Eq. ( 11) denotes thermal diffusivity.The solution of Eq. ( 11) is obtained for the following boundary conditions and the initial condition The problem of ( 11)-( 13) with the initial condition expressed in (14) above is solved using the Laplace transform.The result is ( ) where n = 1, 2, ... and Fo = κt/s 2 designates the Fourier number.
All calculations are performed using the FORTRAN programming language.The calculation results of timedependent temperature changes on both surfaces of the plate obtained using the analytical method are presented in Figure 6(a).Because the plate is heated evenly, the calculated temperature is the same in 25 nodes located on the insulated outer surface.
Then, on the basis of the obtained temperatures on the insulated surface, the temperature inside the wall at z = Δz and z = 2Δz and on the plate inner surface are determined using the inverse marching method.The calculation results are presented in Figure 6(b).The comparison between the temperature values on the plate inner surface obtained by means of the direct method and the inverse marching method is shown in Figure 7. Due to high thermal inertia of the thick-walled plate, the inner surface temperature calculated using the inverse method, with a small error compared to the temperature calculated using the direct method, is obtained with some time lag.However, when a small increase in the outer plate surface temperature is recorded, the inverse method enables rapid determination of the inner surface temperature.Better results can be achieved by carrying out the temperature measurements inside the plate wall, drilling holes in it, but this is not always allowed for components operating under high pressure.

. Conclusions
This paper presents a general method for determination of temperature, the heat flux density and the heat transfer coefficient on the inner surface of a thick-walled element.Based on the described procedure, a calculation program was compiled and validated for determining the plate inner surface temperature.
The inverse method was validated computationally.The comparison between the temperature values obtained from the direct heat conduction problem solution and the results achieved by means of the proposed inverse method is very satisfactory.Stable and accurate measurement results are obtained using the proposed method based on the solution of the inverse heat conduction problem due to partial elimination of random errors in the input signal by using a 9-point moving digital filter [13].
After some modifications, the developed method can be used to monitor thermal stresses and evaluate the degree of wear and estimate residual life of the power plant pressure components such as boilers or turbines.By measuring the pressure element outer surface temperature in a small area, it is possible to determine the temperature of the inner surface.The presented method combined with the method of measuring the fluid unsteady-state temperature described in [12] enables calculation of the heat transfer coefficient on the inner surface.Owing to precise determination of the heat transfer coefficient on the inner surface, thermal stresses arising in the pressure element can be calculated by means of the finite element method.The presented method can be used in online applications.Another advantage of the method is the ease of practical application.

Fig. 1 .
Fig. 1.General view of the plate Based on measured temperature Tm on the plate outer surface (b) shows the cell around the internal node; such a case occurs in the plate distribution into a larger number of cells in the z-axis.All the cells have full dimensions: Δx, Δy, and Δz.The balance equation is written for the central cell.In the case of the heat balance compilation for a node situated at distance Δz from the cooled surface, the top cell has thickness Δz/2, and the rest are full size (Figure 3(c)).

Fig. 3 .
Fig. 3. Diagram illustrating control volumes around the node for which the heat balance equation is written: (a) node located on the plate insulated surface; (b) node located within the space between the plate insulated and cooled surfaces; (c) node located on the cooled surfaceThe heat balance equations are written for all nodes Pi, j, k = P(xi, yj, zk) located in the centres of gravity of each control volume.The coordinates of nodes in the centre of the control volumes are defined by the following formulae: ) and 4).Then time-dependent temperatures in the nodes shown in Figures 3(b) and 4 lying in plane z = Δz are determined from the heat balance equations which are written for the nodes shown in Figures 3(a) and 4. In a similar way, time-dependent temperatures in the nodes in plane z = 2Δz (Figures 3(c) and 4) are determined from the heat balance equations for the nodes lying in plane z = Δz.Next, the temperature in node (4,4,4), lying on surface z = 3Δz = s is determined from the heat balance equations for the nodes lying in plane z = 2Δz.The finite volumes, in the middle of which temperature is measured or determined, lying in different planes z are marked with different colours (Figures 2, 4 and 5).

Fig. 6 .Fig. 7 .
Fig. 6.Time-dependent changes in temperature: (a) data generated by the direct method; (b) results of calculations using the inverse marching method