Inverse Heat Conduction Solution Utilizing the Difference Method with an Exact Matching Rule A Robustly Stable and Easy-to-Use Scheme

Thermal Science & Engineering Vol.29 No.2 (2021). - 33 -. Inverse Heat Conduction Solution Utilizing the Difference Method with an. Exact Matching Rule –– A Robustly Stable and Easy-to-Use Scheme. Yoshihiko HARAMURA 1. Abstract. This paper deals with solutions of inverse heat conduction problems to specify surface heat flux from temperature histories in a body. Heat flux at each timestep is determined as the calculated temperature at a measurement point exactly matching the measurement. This is applied not only to planar bodies, but also cylindrical and spherical bodies. Although this method can cause numerical instability for small timesteps, we find that this is not always true. A fully implicit finite difference method gives stable solutions for any timestep length though there exists restriction of node spacing. Stability of this scheme is explained using the eigenvalue of the matrix. The obtained heat flux waveform shows some delay. We find that a scheme using interpolation with quadratic functions to find the mean temperature of control volumes leads to a delay of at least half a timestep. This interpolation causes numerical instability for very large node spacings, but this restriction is sufficiently acceptable. The effects of noise are studied to show the amplifying factor from temperature to heat flux as a function of timestep. Simulations confirm the effectiveness of this factor. It is also shown that filtering prior to this inverse solution is effective.. Keywords: Heat conduction, Inverse problem, Function specification, Exact matching, Fully implicit. method, Numerical analysis. Nomenclature A : Cross-sectional area [m2] (A–1B)j,j : j-th diagonal element of matrix A–1B b : Thickness or outer radius [m] c : Specific heat [J/(kgK)] C, C0, C1 : Constants EV : Eigenvalue of matrix A–1B [ – ] J : Number of control volumes [ – ] K : Node spacing ratio = ∆zj+1/2/∆zj–1/2 [ – ]. q : Heat flux [W/m2] q1 : Reference heat flux value [W/m2 ] r : Radial coordinate [m] R : Function of r used in separation of variables t : Time [s] T : Function of t used in separation of variables U : Transformed variable = r [Km] Vj : Volume of control volume j [m3]. Y n : Given temperature at adiabatic face [K]. z : Rectangular coordinate [m] Z : Function of z used in separation of variables b j ±  : Weight of neighboring node given by. Eq. (6) (a,b) [ – ]. g j ±  : Parameter given by Eq. (10) [ – ]. ∆r : Node spacing [m] ∆t : Timestep [s]. ∆z : Node spacing [m] ∆j : Dimensionless node spacing = ∆r/rj [ – ]. ∆ : Reference temperature difference = 1 K [K]. ∆ : Amplitude of noise [K].  : Parameter defined with Eq. (32) [– ].  : Time-difference weighting factor [ – ].  : Thermal conductivity [W/(mK)] q j. : Representative value of q j n . between times n and n+1 [K]. q j n : Temperature at location j and time n [K]. : Column vector of q j n [K]. q j. n : Mean temperature giving enthalpy in control volume. [K].  : Density [kg/m3]  : Angular frequency [rad/s]. Suffix j : Space index max : Maximum n : Timestep  : Differential operator + : Dimensionless for t or q. 1 Introduction . It is often necessary to estimate surface heat flux from. Manuscript accepted: September 24, 2020; Editor: Akihiko HORIBE 1 Kanagawa University, Department of Mechanical Engineering (3-27-1 Rokkakubashi, Kanagawa-ku, Yokohama 221-8686, Japan). Thermal Science & Engineering Vol.29 No.2 (2021). - 34 -. internal temperatures. It is widely known that temperature. distributions changing with time can be determined from. the heat conduction equation in a time–space domain after. a time when the temperature distribution is given and in. the space surrounded by a closed surface on which a. boundary condition such as temperature or heat flux is. imposed. (This is called the direct problem, in contrast. with the inverse problem mentioned below.) However, the. problem raised at the beginning is not in this category.. Such problems, called inverse problems, are ill-posed, and. their solutions can be obtained only under some. restrictions.. Inverse heat conduction problems are classified either. as function specifications that find heat flux and/or. temperature histories, or as parameter estimations that. find thermal properties such as thermal diffusivity. This. paper deals with only the former type. . Solving this problem has a rather long history. Stoltz [1]. first proposed a method for obtaining heat conductance on. the surface of a silver sphere quenched in oil. That method. successively finds a new stepwise surface heat flux that. fits the convolution of internal temperature responses. against heat fluxes at instances prior to the measurement.. Stoltz reported that this method led to oscillations for. small timesteps. Beck [2] introduced least-square. regularization to avoid instability. This scheme. determined new surface heat flux, as calculated internal. temperatures match the measurements in an averaged. manner. This regularization also suppresses the effect of. noise superposed on the measurements. Stoltz’s scheme,. in contrast, can be called exact matching.. Many schemes were proposed following those for. estimating surface heat flux [3–13]. Most employed some. form of regularization, making them rather complicated.. Exceptions are Liu et al. [6] and Huang et al. [11], who. adopted exact matching. Both transform the unsteady. conduction equation to two first-order differential. equations introducing V=∂/∂z, which makes the scheme. somewhat complicated. Liu et al. [6] used only an explicit. difference form for time and found instability at small. timesteps identical to a common numerical instability.. Huang et al. [11] introduced a time-difference. weighting factor () and reported that the fully implicit. method has no timestep restriction, the same conclusion. as this paper. However, this does not agree with the. statement that exact matching schemes cause instability. for small timesteps by Beck [2] or Beck el al. [14], which. has long been believed by many researchers, including the. present author. Huang et al. did not discuss estimation. error, either that due to the selection of node spacing or. due to noise. We cannot therefore conduct their scheme. with proper parameters. The objectives of this paper are to. clarify why the implicit method produces stable solutions,. how we should divide the space, and what timestep to. select. We found some differences in solution stability and. in heat flux reconstruction when using a quadratic. interpolation of temperatures or simple node temperatures. to calculate control-volume enthalpy, so we will also. report on these trends.. 2 Formulation. 2.1 Region to be solved and boundary conditions. This study considers only one-dimensional cases,. including planar bodies, long cylinders, and spheres. In. the cases of cylinders or spheres, the temperature sensor. is located at the center. In addition, we assume constant. properties.. Most function specifications, or specifications of. surface heat flux, are performed from temperature. histories at two internal points (letting, say, P1 be the. shallower point and P2 the deeper point; the temperature. at the deeper point may be replaced with a given heat flux. there). Section 1.2 of Beck et al. [14] points out that this. situation is divided into two regions, treated as direct and. inverse problems, because the temperature distribution. with time in the region between P1 and P2 can be solved,. and the heat flux at P1 is also obtained. Thus, the region. to be solved as a pure inverse problem is the region. shallower than point P1. Temperature and heat flux are. both given at one boundary z= 0, and the surface is located. at z=b as shown in Fig. 1. Heat flux at the surface will be. estimated.. 2.2 Energy equation and its boundary conditions. We apply the finite difference method. Nodes are. Fig. 2 Node arrangements.. Fig. 1 Region to be solved for the constant cross-. section case.. Thermal Science & Engineering Vol.29 No.2 (2021). - 35 -. located at the circles in Fig. 2. The first node (j = 1) is located at the measurement point z = 0 or r = 0. The last node (j = J+1) is a virtual node located outside of the body, as the surface is at the midpoint between this node and its neighboring interior node (j = J). Control volume j is bounded by planes that cross the midpoints of node j in the volume and the neighboring nodes, j–1 and j+1. This arrangement gives the highest accuracy when describing heat flow across the boundary.. The basic energy equation is . rcV j. q j n+1. -q j n. ∆t = A. j-1 2. l q j-1. -q j. ∆z j-1 2. + A j+1 2. l q j+1. -q j. ∆z j+1 2. (1a). for control volume j (from 2 to J) The superscript denotes the timestep and the subscript is the index of the node or the control volume. q j. n is the mean temperature in. control volume j. q j without an index n denotes some representative value of the node temperature at a time between times n and n+1. In the special case of a planar body control volume in 0 ≤ z ≤ ∆z1+1/2/2, the equation corresponding to Eq. (1a) is . rcA ∆z. 1. 2. q 1. n+1 -q. 1. n. ∆t = Al. q 2 -q. 1. ∆z 1+1 2. , (1b). because the left boundary is adiabatic. This condition is due to the test case explained below. For a control volume in 0 ≤ r ≤ ∆r/2,. rcp ( ∆r 2. ) 2 q. 1. n+1 -q. 1. n. ∆t = 2p ( ∆r. 2 ) l. q 2 -q. 1. ∆r (1c) . is satisfied for a unit length of cylinders and. rc 4p. 3 ( ∆r 2. ) 3 q. 1. n+1 -q. 1. n. ∆t = 4p ( ∆r. 2 ) 2. l q 2 -q. 1. ∆r (1d) . for spheres. Here . ∆z j-1 2. = z j - z. j-1 , ∆z. j+1 2 = z. j+1 - z. j , (2). Aj±1/2 are the cross-sectional areas at the boundaries of a. control-volume j. They are constant A for planar bodies. The volume Vj is described as Aj ∆zj for planar or. cylindrical bodies, with . ∆z j = ì í ï. îï. 1. 2 ( z j+1 - z j-1 ) (for j ³ 2). z j+1. - z j. (for j =1). (3). and the cross-sectional area Aj at a node j, and. V j =. 4p. 3 ( r. j+1 2. 3 - r. j-1 2. 3 ) (4). for spheres. In the case of a cylinder or a sphere, z in Eqs. (1a), (2) and (3) is replaced with r. q j. n , the mean. temperature giving enthalpy in the control volume, is. q j. n = b. j -q. j-1. n + (1- b. j - - b. j + )q. j. n + b. j +q. j+1. n , (5). where b j ±. for variably spaced nodes in planar bodies are. b j ± =. (∆z j±1 2. )2 + 2∆z j+1 2. ∆z j-1 2. - 2(∆z j 1 2. )2. 24∆z j ∆z. j±1 2. (6a). for j ≥ 2 or. b 1. +. =1 12 (6b). for j = 1. Table 1 lists values b j ±. for cylinders or spheres with equally spaced nodes. b1. -. does not exist. The quadratic temperature distribution in each control volume against the location is postulated in the derivation of. b j ± . including j = 1, where the gradient is zero at z = 0 or r = 0. q j in the conduction terms in Eqs. (1a)–(1d) is expressed using a time-difference weighting factor  as . q j = hq. j n+1. + (1 -h)q j n . (7). Parameter  is 0.5 for the Crank–Nicolson scheme and 1 for the fully implicit scheme.. Equations (1a)–(1d) thus become. -C j - q j-1 n+1. + (1+C j - +C. j + )q. j n+1. -C j + q j+1 n+1. = D j - q j-1 n. + (1-D j - -D. j + )q. j n +D. j + q j+1 n. (8). with. C j ± = hg. j ± - b. j ± , D. j ± = (1-h)g. j ± +b. j ± (9). and. g j ± =. l ∆t. rc∆z j ∆z. j±1 2. . (10). Table 1 also lists values g j ±. for cylinders and spheres with equally spaced nodes. We often use. q j. n = q. j. n (11). instead of Eq. (5) with Eqs. (6)(a,b) for simplicity, in which case b j. - and b j. + will be zero. We call the scheme. using Eq. (5) with Eqs. (6)(a,b) Scheme A and that using Eq. (11) Scheme B. Unless otherwise noted, we employ Scheme A.. Table 1 Parameters in Eq. (9) for equally spaced nodes. Cylinder Sphere. b j ±. (j≠1) 1. 24 (1 ± ∆x j ). 1. 24 ± 1. 12 ∆x. j +. 1. 160 (∆x. j )2. b 1. + 1. 8. 3. 20. g j ±. (j≠1). l ∆t. r c(∆r)2 (1 ± 1. 2 ∆x. j ) l∆t. r c (∆r)2 (1 ± 1. 2 ∆x. j ) 2 . g 1. + 4 l∆t. r c(∆r)2 6. l ∆t. r c(∆r)2. ∆x j = ∆r r. j. ±. Thermal Science & Engineering Vol.29 No.2 (2021). - 36 -. The heat flux boundary condition is already given as. Eqs. (1b)–(1d). The remaining boundary condition is that. for temperature at z= 0 or r= 0, where temperature q 1. n+1. equals the given temperature Yn+1, that is,. q 1. n+1 =Y. n+1 . (12). This is the condition of exact matching. . 2.3 Coefficient matrices of linear system for q j n+1. . Equations (8) and (12) form a linear system of J+1. algebraic equations with J+1 unknowns,. . (13). A, B, and b (a column vector) can be arranged as. (14). (15). and. . (16) . The first row represents boundary condition (12).. Multiplying A–1 from the left, Eq. (13) becomes. . (17). This equation is repeatedly applied to evolve the. temperature distribution. The surface heat flux added at. time n+1 is calculated as q = l (q J+1. n+1 -q. J. n+1 ) ∆z. J+1 2 .. 2.4 Given temperature on adiabatic face or origin. Temperature response at the adiabatic face z=0 or the. center r= 0 against the prescribed heat flux is used as the. given temperature. We select the heat flux waveform. following the test cases of triangular heat flux presented. by Beck et al. [14]. The initial condition is q j 0= 0 for all. locations. Figure 3 shows the heat flux history. Heat flux. increases linearly until dimensionless time t+, defined as . t +. =. l. rcb2 t , (18). reaches 0.6, and then linearly decreases to zero by t+= 1.2.. Because these periods are too short for the large ∆t+. recommended for use with the exact matching scheme, we. also use 1.2 and 2.4 in place of 0.6 and 1.2 in sections 3.3. and 3.4 for ∆t+> 0.07. Dimensionless heat flux is defined. as. q + =. q. q 1. , (19). where. q 1 =. l. b ∆q . (20) . Here ∆ is the reference temperature difference,. considered to be 1 K. . The given temperatures Y n are calculated as those for. the direct problem, solved by the difference method of the. Crank–Nicolson scheme with J = 29 and ∆t+= 1.010–4.. The error is estimated as less than 0.01% of ∆, from. comparison with an analytical solution [15] for a constant. heat transfer coefficient l b (Biot number = 1) boundary. condition instead of the heat flux input described in Fig. 3.. Analytical (z = 0) [14] was also used for planar bodies,. and the results are the same, within 0.1% of the maximum. value of q+.. The calculation was performed using dimensional. quantities. Settings b = 5 mm, = 400 W/(mK),c. = 4.0106 J/(m3K) are used for both the direct and inverse. solutions. . Fig. 3 Heat flux history of test case.. 3 Results and discussion. 3.1 Stability of solutions for planar bodies. This section considers problems in which surface heat. flux is estimated from the temperature at an adiabatic face. of planar bodies. Figure 4 shows results obtained for. ∆t+= 0.01. The solution is quite stable when  is near unity. (over 0.9 in this case), though there is some delay and. rapidly decaying fluctuation can be observed just after. timings where the original heat flux varies its changing. rate (t+=0, 0.6 or 1.2).. Beck el al. [14] showed unstable results for ∆t+= 0.05. (stable for ∆t+≥ 0.3) calculated through the Stoltz (exact. matching) method. This is also true under the Crank–. Nicholson scheme (= 0.5) in the present calculation.. However, most results using the fully implicit scheme are. Thermal Science & Engineering Vol.29 No.2 (2021). - 37 -. stable. This stability can be theoretically explained as follows.. Coefficient matrices A and B are lower triangular ones, so A–1B also becomes a lower triangular matrix, the diagonal elements of which are. (A-1B) j+1, j+1. = -. 1-h + b j + g. j +. h - b j + g. j +. (21). except for the first, corresponding to the temperature. boundary conditions, which is zero. Here b j + g. j +. is. b j +. g j +. =. 1. 24 ( ∆z. j+1 2. b ) 2 1. ∆t +. (22a). for equally spaced nodes or. b j +. g j +. =. K 2 + 2K - 2. 24K 2. ( ∆z. j+1 2. b ) 2 1. ∆t +. (22b). for variably spaced nodes, where K is the node spacing. ratio ∆zj+1/2/∆zj–1/2. The diagonal elements of a triangular . matrix are also eigenvalues of the matrix. As the evolution. proceeds, each element of an eigenvector or each mode. will change as its eigenvalue is multiplied. Since the. actual temperature distribution is considered to be the sum. of the modes, the existence of at least one eigenvector that. diverges or converges too slowly may spoil the solution.. When ∆t+ is large and/or ∆z is small as b j +. g j +. . becomes negligible, eigenvalues become -(1-h ) h . If.  = 1 (the fully implicit scheme), eigenvalues are all zero,. so any disturbance will immediately disappear. If,. however, = 0.5 (the Crank–Nicolson scheme),. -(1-h ) h becomes –1 and the disturbance persists in. an oscillating manner, as its amplitude remains constant.. Haramura [16] showed that the amplitude does not go to. infinity, instead approaching a large constant. Inherent. unstable characteristics of the inverse problem seem. suppressed by the damping property of the fully implicit. scheme. When  has an intermediate value, the. disturbance will decay at some timesteps. Figure 4 shows. this damped oscillation for the case of = 0.75 as a thin. broken line. There is no spatial wave even in such. oscillating periods.. Next, we discuss the stability for a nonnegligible. b j +. g j +. . When. b j +. g j +. is greater than 1/2, eigenvalues. exceed –1 even for = 1, so the solution becomes unstable.. A smaller b j +. g j +. is preferable for fast damping. When. the maximum absolute value of (negative) eigenvalues. EV are prescribed as disturbance decays for the. designated rate, the relation between ∆zj+1/2 and ∆t + is. calculated from Eqs. (21) and (22b) as. ∆z j+1 2. b = é ëê C ( h - 1. 1- EV )∆t+ ùûú. 1 2. , (23). where. C = 24K. 2. K 2 + 2K - 2. . (24). When j =1, C will be 24 for both equally and variably. spaced nodes. Equation (23) gives the maximum node. spacing for a stable solution. Since the maximum ∆z. causes instability, equally spaced nodes are best for the. smallest number of divisions. When = 1, EV = –1 3, and. K = 1, the factor in front of ∆t+ becomes 6. ∆z/b less than. 0.24, that is, J ≥ 5 (b ∆r = 4.5 when J = 5), will stabilize. the solution even when ∆t+ is as small as 0.01.. To confirm the stability criteria, Fig. 5 shows ∆z/b and. ∆t+ combinations for stable (circles) and unstable. (crosses) solutions for equally spaced nodes. Here,. stability was judged using an amplifying factor by which. the peak-to-peak heat flux value grows at each timestep.. This peak-to-peak value is calculated from the difference. in heat flux at a given time and the average of those one. step earlier and later, as shown at the right lower corner of. the figure. The period for finding changes in peak-to-peak. values was ten timesteps just after t+= 1.7. When variably. spaced nodes are used, we can easily find the neutral. Fig. 4 Effect of time-difference weighting factor  for. planar bodies. . Fig. 5 Upper limits on node spacing giving stable. solutions for planar bodies. The figure at the lower right. shows how to find peak-to-peak values to calculate. amplifying factors.. Thermal Science & Engineering Vol.29 No.2 (2021). - 38 -. stable condition at which the peak-to-peak value of heat. flux fluctuation remains constant. Figure 5 also shows. combinations of ∆zmax/b (∆zmax is the largest ∆zj+1/2) and. ∆t+ of the neutral stable state, with squares for K<1 at. which the node spacing ∆z1+1/2 at the adiabatic face. controls the stability, and marks. ,. , and . for K>1 at. which the node spacing ∆zJ+1/2 at the surface controls. stability. Geometric series of node spacings are employed. here.. Circles are located below and crosses are above Eq.. (23), with EV = –1 and K= 1 plotted as a thick solid line.. Squares are located on that line. The other marks for K>1. are located on different lines, depending on the number of. divisions J, which determines different K for a given. ∆zJ+1/2 because of the restriction . . (25). All are located above Eq. (23) with EV = –1 and K= 2, at. which C takes a minimum of 16. These results support the. validity of criteria, Eq. (23). . There are two other possibilities for this formulation.. One is to use Eq. (11), Scheme B, instead of Eq. (5) with. Eqs. (6)(a,b). This makes b j +. g j +. zero, then eigenvalues. -(1 -h ) /h . This solution is always stable, independent. of node spacing, but section 3.3 will show one demerit of. Scheme B. Another possibility is to arrange node locations. so that the adiabatic face is located at the midpoint. between the two leftmost nodes, the same as the surface.. This eliminates the different formula of Eqs. (1b) and (6b). for the first control volume, slightly improving stability. but having the same demerit as Scheme B if a quadratic. interpolation of temperature is employed. . Summing up, the best selection for a stable solution is. to use equally spaced nodes satisfying Eq. (23) with C =24. and EV ≥ –1 3 or so. Any merits of using  other than. unity seem negligible.. 3.2 Results with cylindrical or spherical bodies. This section considers problems in which surface heat. flux is estimated from the temperature at the center of a. cylinder or a sphere. Figure 6 shows some examples of. heat flux waveforms with a result from a planar body.. They are quite identical. Stability can be similarly discussed. Matrix A–1B will. be a triangular matrix, with eigenvalues again given by Eq.. (21). Since equally spaced nodes are employed, ∆x j in. Table 1 is 1 ( j -1) , and b j +. g j +. is calculated as. b j +. g j +. =. 1. 12. j. 2 j -1 ( ∆r b. ) 2 1. ∆t+ (j ≥ 2), (26a). except for j = 1, where the factor in front of (∆r b)2 ∆t+ . is 1 32 for cylinders, or as. b j +. g j +. =. 1. 6. j2 -17 20. (2 j -1)2 ( ∆r b. ) 2 1. ∆t+ (26b). including j= 1 for spheres. From Eqs. (26)(a,b) including. j = 1, b 2. +. g 2. +. is largest as . ( b j. +. g j. + ) max. =. b 2. +. g 2. + =. 1. 18 ( ∆r b. ) 2 1. ∆t +. (27a). for cylinders and. ( b j. +. g j. + ) max. =. b 2. +. g 2. + =. 7. 120 ( ∆r b. ) 2 1. ∆t +. (27b). for spheres. From Eq. (21) and Eq. (27a) or (27b), ∆r for. which all eigenvalues are smaller (nearer 0) than a. prescribed negative eigenvalue EV is given by Eq. (23). with C = 18 for cylinders and C = 120/7 for spheres. Here. ∆zj+1/2 in Eq. (23) are replaced with ∆r.. Figure 7 shows the solution stability as judged from the. amplifying factor of heat flux fluctuation mentioned in. section 3.1, with diamonds and circles respectively. showing combinations of ∆r/b and ∆t+ for which the. solution is stable (the amplitude of fluctuation decays) for. cylinders and spheres. Crosses indicate ∆t+ at which the. Fig. 6 Comparison of heat flux waveforms for. cylindrical, spherical, and planar bodies.. Fig. 7 Upper limit on node spacing for cylinders and. spheres. Numeric labels are amplifying factors of heat. flux fluctuation.. Thermal Science & Engineering Vol.29 No.2 (2021). - 39 -. solution loses stability for the same ∆r. Numeric labels are. amplifying factors at each timestep. Factors especially. near criteria are almost the same as eigenvalues.. Underlined factors are not, but fluctuations for these. conditions are already damped to less than 110–5 in the. dimensionless amplitude, so factor accuracy is doubtful.. We can therefore conclude that. b 2. +. g 2. +. , the absolute. maximum of eigenvalues, controls the stability.. 3.3 Delay. As mentioned in the last part of section 3.1, the solution. using Eq. (11), Scheme B, is stable even for large ∆z. We. therefore first examined reconstruction of heat flux with. this scheme, using different values for the numbers of. divisions J in a range from 19 to 2 using equally spaced. nodes. The relation between J and ∆z b is . ∆z. b =. 1. J -1 2 . (28). Figure 8 shows some results, which are nearly the same. except for the delay. The plots for J =5 and J =19 are. hardly distinguishable; only a very small J makes the. delay large. . The delay value is estimated by fitting a linear function. with the same gradient as the original in the region. 0.2 ≤ t+≤ 0.6 (0.4 ≤ t+≤ 1.2 for the slow heat flux pattern. employed for ∆t+> 0.07). These are plotted in Fig. 9, with. solid symbols denoting Scheme A in which Eq. (5) with. Eqs. (6)(a,b) was used, and open symbols denoting. Scheme B in which Eq. (11) was used. . Delay is large for large ∆t+. As J decreases, the delay is. constant in the case of Scheme A, but significantly. increases under Scheme B over Scheme A values. Despite. the smaller delay, Scheme A tends to become unstable for. some large node spacings and short timesteps, as shown. in section 3.1. . Note that the delay under Scheme A and the lower limit. under Scheme B is near ∆t+ 2 . This implies that the. fully implicit scheme will lead to a half timestep delay due. to its difference method, while the central difference. method with = 0.5 leads to zero delay. That the delay for. = 0.75 is nearly half of that for  = 1 supports this. interpretation. . Since calculation efficiency (number of arithmetic. operations) does not improve, there is no merit to. employing Scheme B, including the cases of cylinders or. spheres, where delay remains as ∆t+ 2 so long as we. use. b j ± and. g j ± in Table 1. . 3.4 Effect of noise. The most severe restriction on heat flux specification is. caused by noise superposed on the temperature. measurement. This depends on the measurement. technique, so we focused on sensitivity factors from . temperature to surface heat flux q q r or z=0. . These are. obtained from the frequency response for Nyquist. frequency 1 (2∆t), expressed in complex form as . b. l. q. q z=0. = (1+ i)e ( sinhe cose + icoshe sine ) (29). for planar bodies, . b. l. q. q r=0. = 2e ( be ¢r 2e + i be ¢i 2e ) (30). for cylinders, and. b. l. q. q r=0. = 1. 2. é ëê ( coshe -. 1. e sinhe )cose. + ( cose - 1 e sine )coshe + i{ 1. e ( sinhe cose. - coshe sine )+ 2sinhe sine } ùûú. (31). for spheres, where. Fig. 9 Delay of heat flux in the first half of the heating. period.. Fig. 8 Effect of number of divisions J on heat flux. waveform for planar bodies with equally spaced nodes.. Thermal Science & Engineering Vol.29 No.2 (2021). - 40 -. e = rcb2w. 2l =. p rcb2. 2l∆t . (32). ber(x) and bei(x) are Kelvin functions [17] and prime (). denotes derivative d/dx. The appendix shows derivations. of Eqs. (29) – (31).. Figure 10 shows absolute values of Eqs. (29) – (31). For. all shapes, the factor is inversely proportional to ∆t+ for. large ∆t+>1. One half or one third of ∆t in planes is. available for cylinders and spheres in this range, as can be. seen from the formula for the asymptotes indicated in the. figure. For smaller ∆t+, the factors are larger than the. asymptotes.. For copper bodies (assuming = 400 W/(mK), c = 4.0. 106 J/(m3K)) with b = 5 mm, noise of 0.01 K at the. Nyquist frequency will be amplified to 6.1, 2.7, or 1.8. kW/m2 when ∆t+= 0.5, and 26, 9.8 or 5.6 kW/m2 when. ∆t+= 0.2 for planar, cylindrical, or spherical bodies,. respectively. Although the limitation of ∆t+ depends on the. actual noise level, ∆t+= 0.2 or less may be applicable even. in the exact matching scheme.. Figure 11 shows examples of heat flux waveforms. affected by noise. Noise superposed on the given. temperature is regular noise, which alternates every. timestep with a constant amplitude, except that it is zero. at t+= 0. Figure 10 shows four ∆t+ conditions (∆t+= 0.0887. and 0.1612 for a plane, 0.0977 for a cylinder and 0.0708. for a sphere), using the same symbols as in Fig. 11.. (b l )q q r or z=0. = 200 for the first condition and 50 for. the other three. The amplitudes of noise ∆N are selected. so that the amplitude of q+ will be 0.2 for all conditions if. Eqs. (29) – (31) work. This amplitude is indicated as two. thin broken lines. The actual amplitude is around 1.3 times. this theoretical value. This discrepancy is probably due to. a small transient overshoot, as observed in Figs. 4, 6, and 8. . Most practical solutions of inverse heat conduction. problems adopt regularization to reduce the effect of noise,. so it is ineffective to apply additional filtering before such. process. Now, however, we have acquired a stable, exact. matching scheme, so applying a low-pass filter prior to. processing seems effective. The thick broken line with. circles in Fig. 11 indicates the heat flux obtained through. such processing. The given temperature superposed with. noise is the same as the case with squares. A Savitzky–. Golay filter [18] with M = 2 (to approximate local. quadratic functions), nL= 2, and nR= 2 (using data from. 2nd earlier to 2nd later) was applied before solving. This. filtering reduces the effect of noise to about one third. without significant distortion. More effective filters. should result in further reduction.. Fig. 11 Response to temperatures with regular noise. for planar (��), cylindrical (�) and spherical (�). bodies. � used filtered data of �.. 4. Determining timestep and node spacing. 4.1 Solid cylinders or spheres. To estimate surface heat flux from a measured. temperature at the center of a solid cylinder or sphere (see. section 4.2 when dealing with hollow cylinders or. spheres), equally spaced nodes are preferable for. simplicity, since the eigenvalues of A–1B are somewhat. similar. We cannot specify the exact initial condition but. the difference from it exponentially decays over time. The. solution to the Capstone problem [15] shows that such a. difference decays to less than 1% of the initial value at. t+= 2, assuming thickness 2b. It is sufficient to use the. steady solution that satisfies the two measured values at. time zero or the linear distribution that passes the two. measurement values at time zero, . The procedure for finding a temperature distribution. that satisfies the boundary conditions at two locations at. every timestep is as follows.. (1) Select a noise reduction filter. A noise reduction filter is essential to make. timestep ∆t small enough to cope with quickly. changing phenomena. Powerline frequency and the. Nyquist frequency should be focused on. . (2) Determine timestep ∆t. ∆t+ then ∆t are determined from Eq. (30) or (31). from the noise level of measured temperature after. Fig. 10 Sensitivity factor from temperature at the. adiabatic face to surface heat flux.. Thermal Science & Engineering Vol.29 No.2 (2021). - 41 -. filtering and the required certainty of heat flux. If ∆t. is prescribed from the experiment, an appropriate. filter should be selected to satisfy the required. certainty of heat flux.. (3) Determine node spacing ∆r. Setting = 1 and EV = –1 3 or so in Eq. (23) with. C = 18 or 120/7, find the maximum ∆r, then determine. the number of divisions J. Follow the node. arrangement rule presented in Fig. 2 for higher. accuracy.. 4.2 Planar bodies. When we estimate surface heat flux on planar bodies,. there are three characteristic locations: the surface, the. shallower temperature sensor, and the deeper sensor. It is. best for accuracy to select node location such that the. measurement points match the nodes, and a virtual outside. node set as the surface is at the midpoint between this and. the neighboring inside node, as in Case I in Fig. 12. The. third point, the location of deeper sensor, may be replaced. with a face where known heat flux is applied. In that case,. the node arrangement should be changed as in Case II in. Fig. 12. . In either case, node temperatures in the direct region,. the region between two measurement points, are. determined independently of temperatures in the inverse. region, because the number of control volumes plus two. conditions (corresponding to boundary conditions) are. given for the same number of unknowns. As for node. temperatures in the inverse region, A–1B is a lower. triangular matrix of diagonal elements given with Eq. (21).. The stability criteria are therefore given by Eq. (23) with. C = 24 for a designated EV. See the previous section. regarding the initial condition.. Follow the procedure from (1) to (2) in the previous. section and (3) below.. (3) Determine node spacing ∆z. Set the node spacing in the inverse region uniform to. make all eigenvalues in the inverse region the same. This. results in the fastest decay of disturbance for a fixed J.. Find ∆zj+1/2 in the inverse region to satisfy the stable. condition of Eq. (23) with C = 24,  = 1 and EV = –1 3. or so. The same ∆zj+1/2 is applied to the node next to the. shallower measurement point in the direct region. . There is no restriction on node spacing in the direct. region, so it is sufficient to select the different uniform. spacing in the direct region with a spacing near ∆zj+1/2, as. above. . 5 Conclusion. We applied a finite difference method to an inverse heat. conduction solution in which the temperature at. measurement point exactly matches the measurement.. Primary findings were as follows:. (1) A fully implicit method gives a stable solution for any. timestep, even in the case of exact matching. This is. valid not only for planar bodies, but also for. cylindrical and spherical bodies.. (2) Estimated heat flux is delayed by at least ∆t 2 from. the original. This is probably due to inherent. characteristics of the fully implicit method. This. minimum delay is maintained for wider node spacing. when quadratic interpolation is used to find the mean. temperature in control volumes. The delay increases. with wider node spacings when the node temperature. is considered as the mean temperature.. (3) The scheme using quadratic interpolation of. temperatures becomes unstable when the node. spacing becomes large. The optimum node spacing in. the pure inverse region is given by Eq. (23). This. should be used with = 1 and EV  –1 3.. (4) We estimated the effect of noise with the frequency. response for the Nyquist frequency. Actual amplitude. is about 1.3 times the theoretical value. The exact. matching scheme can be applicable at ∆t+= 0.2 or less.. We showed that low-pass filtering prior to the. solution is effective.. (5) We presented guidelines for determining timesteps. and node arrangements for actual cases where surface. heat flux is estimated from two measurements in a. body.. We acquired a robustly stable scheme with the. negligible sacrifice of a ∆t 2 delay caused by fully. implicit differencing. This scheme is computationally. quite effective and easily implemented. It can also be. easily extended to multidimensional problems to evaluate. changes in surface heat flux distributions over time, when. the side faces can be treated as adiabatic. It remains as a. future task to evaluate the effect of deviations of assumed. temperatures or heat fluxes on the side faces from the. exact distributions.. Acknowledgement. We thank Prof. M. Homma of Kanagawa University for. the advice on matrix eigenvalues.. Reference. [1] Stoltz, G., Jr., “Numerical solution to an inverse. 1. Fig. 12 Node arrangement for actual cases.. Thermal Science & Engineering Vol.29 No.2 (2021). - 42 -. problem of heat conduction for simple shapes”, J. Heat Transfer, 82 (1960), 20-26.. [2] Beck, J. V., “Surface heat flux determination using an integral method,” Nuclear Engineering and Design, 7, (1968), 170-178.. [3] Monde, M., “Analytical method in inverse heat transfer problem using Laplace transform technique,” Int. J. Heat Mass Transfer, 43, (2000), 3965-3975.. [4] Woodfield, P. L., Monde, M. and Mitsutake, Y., “Improved analytical solution for inverse heat conduction problems on thermally thick and semi- infinite solids,” Int. J. Heat Mass Transfer, 49, (2006), 2864-2876.. [5] Liu, C-S., “Group preserving scheme for backward heat conduction problems,” Int. J. Heat Mass Transfer, 47, (2004), 2567-2576.. [6] Chang, C-W., Liu, C-S. and Chang, J-R., “A group preserving scheme for inverse heat conduction problems,” Comput. Model Eng. Sci., 10, (2005), 13-38.. [7] Liu, C-S. and Chang, C-W., “A global boundary integral equation method for recovering space-time dependent heat source,” Int. J. Heat Mass Transfer, 92, (2016), 1034-1040.. [8] Khajehpour, S., Hematiyan, M. R. and Marin, L., “A domain decomposition method for the stable analysis of inverse nonlinear transient heat conduction problems,” Int. J. Heat Mass Transfer, 58, (2013), 125-134.. [9] Woodbury, K. A. and Beck, J. V., “Estimation metrics and optimal regularization in a Tikhonov digital filter for the inverse heat conduction problem,” Int. J. Heat Mass Transfer, 62, (2013), 31-39.. [10] Najafi, H., Woodbury, K. A. and Beck, J. V., “Real time solution for inverse heat conduction problems in a two-dimensional plate with multiple heat fluxes at the surface,” Int. J. Heat Mass Transfer, 91, (2015), 1148-1156.. [11] Huang, M. et al., “A new efficient and accurate procedure for solving heat conduction problems,” Int. J. Heat Mass Transfer, 111, (2017), 508-519. [12] Wen, S. et al., “Solution of inverse radiation- conduction problems using a Kalman filter coupled with the recursive least-square estimator,” Int. J. Heat Mass Transfer, 111, (2017) 582-592.. [13] Wen, S., et al., “Efficient and robust prediction of internal temperature distribution and boundary heat flux in participating media by using Kalman smoothing technique,” Int. J. Heat Mass Transfer, 147, (2020), 118851. . [14] Beck, J. V., Blackwell, B. and Clair , C. R. Jr., Inverse Heat Conduction, Ill-posed Problems, (1985), 3-7, 87-94, 169-170, John Wiley & Sons.. [15] Hahn, D. W. and Özisik, M. N., Heat Conduction, 3rd Ed., (2012), 193-196, 208-211, 326-332, John Wiley & Sons, Inc.. [16] Haramura, Y., “Effect of differencing scheme on the stability of inverse heat conduction solution” (in Japanese), Proc. 55th National Heat Transfer. Symposium of Japan, (2018), E123. [17] Moriguchi, S., Utagawa, S. and Hitotsumatsu, S.,. Mathematical Formula III (in Japanese), (1975), 174-177, Iwanami Shoten.. [18] Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P., Numerical Recipes in Fortran 77, 2nd Ed. (1992), 644-649, Cambridge Univ. Press.. Appendix: Derivation of Equations (29)–(31). A solution to a one-dimensional heat conduction. equation in Cartesian coordinates is assumed to be in the form of the product of T(t), a function of time only, and Z(z), a function of position only. When considering a periodic change with an angular frequency , Z(z) must satisfy. ¢T. T =. l. rc. ¢¢Z. Z = iw , (33). so Z is expressed with constants C0 and C1 as. Z =C 0 e (1+i)e z b. +C 1 e -(1+i)e z b. , (34). with. e = rcw. 2l b . (35). Applying boundary conditions . ¢Z (0) = 0 (36). leiw t ¢Z (b) = q eiw t (37). yields. q = 2C 0 e iw t é. ëê coshe. z. b cose. z. b + isinh e. z. b sine. z. b. ù ûú ,. (38). with C0 satisfying. qb. 2lC 0. = (1 + i)eiw t e ( sinhe cose + icoshe sin e ) .. (39) . Dividing Eq. (39) by Eq. (38) and setting z = 0 yields Eq. (29). . For cylindrical coordinates, setting  =T(t) R(r), . ¢T. T =. l. rc ( ¢¢R R. + 1. r. ¢R. R ) = iw (40). is obtained instead of Eq. (33). The function that R(r) satisfies is Kelvin’s [17], and the solution that is finite at r = 0 will be. q =C 0 e iw t é. ëê ber ( 2e r. b ) + ibei ( 2 e r. b ) ùûú (41). with a constant C0 and the same  described in Eq. (35). Combined with the heat flux expression. Thermal Science & Engineering Vol.29 No.2 (2021). - 43 -. q = lC 0 eiw t. 2e. b ( be ¢r 2e + ibe ¢i 2 e ) (42). added on the surface r=b, Eq. (30) is obtained. In the case of spherical coordinates, employing the. variable transformation . U = rq (43). introduced by Hahn and Özisik [15], the heat conduction equation becomes . ¶U. ¶t =. l. rc. ¶ 2 U. ¶r 2. . (44). The solution finite at r = 0 then becomes. q =C 0. e iw t. r. é ëê sinh( e r. b ) cos (e r. b ). + i cosh( e r b ) sin (e r. b ) ùûú. (45). with a constant C0 and the same  described in Eq. (35). This leads to Eq. (31) through a similar arrangement as when going from Eqs. (41) and (42) to Eq. (30).