l.7r.'. . .., .v.-.-.r y \ . .1'1 MS .13.”, _ Ei‘m‘jt ‘ ”if fiwsg' I '1': 34:” 3.3,}; ’"u 255‘“ ”9&3 » 4"". H 42",}: $31)}? v ' 1'}, vii-"u ”on u 1 ..; w.}~:~,:;.».~;-z- ’. :17?) (‘3' 6:21 I 576;: (’fY‘J‘ all} 2...: .._u ./.4’.-.'.«"- '« V. _. . “‘x'v‘fiu‘ ’2‘ ; . «m If ‘3 I] . 4r . ‘-_ . 3-. .u‘.‘\ ‘ {1‘} .IIA' " " , , '. 3‘ .‘ «w :- 7.. -\ta-“ . ‘fi...’ , .uvx {a}? ,3 q?! ‘4‘. -. ._ . -. ._ l 41- l y ‘ F.- maxi “a . ,. ‘~‘ b ‘ a”;- J - v a: n'v‘; Mb 55 NIVEPSIYY LlBRARIES wm‘HMI'imH.Immwnu lfl :H P3 1293 0060000 000 I gdfiESbQ LIBRARY W Michigan State J“Luniversity This is to certify that the dissertation entitled 3o1ution of General Two-Dimensional Inverse Heat Conduction Problems and One-Dimensional Inverse Melting Problems presented by John Siuming Tu has been accepted towards fulfillment of the requirements for Ph.D. degree in Mechanical Engineering gymmi/ Major pfifessor MS U it an Affirmative Action/Equal Opportunity Institution 0-12771 A I: PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE ll ‘ '0” l s. l J MSU Is An Affirmative Action/Equal Opportunity Inmhution SOLUTION OF GENERAL m-DIHENSIORAL DIVERSE HEAT CONDUCTION PROBLEMS AND ONE-DIMENSIONAL INVERSE HELIING PROBLEMS BY John Siuning Tu A DISSERTATION Submitted to Michigan State university in partial fulfillment of the requirements for the degree of DOCTOR OF PHIIDSOPHY Department of Mechanical Engineering 1988 74059 ABSTRACT SOEUTION OF GENERAL TWO-DIMENSIONAL INVERSE HEAT CONDUCTION PROBLEMS AND ONE-DIMENSIONAL INVERSE MEETING PROBLEMS BY John Siuning Tu The first task of this study is to present the solution method for general two-dimensional inverse heat conduction problems. The combined method of function specification and regularization is employed to solve these problems. The piecewise-polynomial function, which can be the piecewise combination of different degree polynomials, is used to approximate the spatial distribution of the unknown heat flux. Any discontinuity in the surface heat flux can also be treated by this piecewise-polynomial approximation. A special programming concept has been employed to develop the inverse problem computer program, IHCPZD, to solve general two- dimensional inverse problems. The strategy is to use an existing code as the direct problem solver. This gives greater power, generality, and reliability for the same effort than if a direct problem code were written specifically for inverse problems. Program IHCPZD has been tested successfully for several cases, and one of these cases is presented. The second task of this study is to solve the one-dimensional inverse melting problem. The interface tracing method is employed to solve the direct one-dimensional melting problem, and program IHCPZD is modified by replacing the direct problem solver to solve this inverse melting problem. Different types of measurements (temperatures and interface location) have been used to estimate the unknown surface heat flux in this problem. ACKNOWLEDGNENTS The author would like to express his appreciation to Dr. James V. Beck for his guidance, support, and encourage on this work. The precious suggestions from Dr. Arafa Osman and the assistance in editing by Ms. Laura Wheeler are also greatly appreciated. In addition, the author would like to thank General Dynamics Convair Division and Acurex Corporation for their financial support on this research. TABLE OF CONTENTS LIST OF TABLES .............................................. vi LIST OF FIGURES ............................................. Vii NOMENCLATURE ................................................ ix CHAPTER 1 INTRODUCTION ...................................... 1 1.1 APPLICATIONS ....................................... 2 1.2 LITERATURE REVIEW .................................. 3 1.3 OUTLINE .......................... i .................. 5 CHAPTER 2 GENERAL DISCUSSION OF THE FUNCTION SPECIFICATION METHOD .................................................. 8 2.1 INTRODUCTION ....................................... 8 2.2 PROBLEM DESCRIPTION ................................ 11 2.3 OUTLINE OF THE FUNCTION SPECIFICATION METHOD ....... 13 2.4 ASSUMING A FUNCTIONAL FORM FOR THE UNKNOWN HEAT FLUX 15 2.4.1 GLOBAL SPECIFICATION VERSUS SUBDOMAIN SPECIFICATION ............................... 1 ...... 15 2.4.2 FUNCTION TYPES .............................. 17 WHOLE DOMAIN ESTIMATION VERSUS SEQUENTIAL ESTIMATION 18 DEFINING THE OBJECTIVE FUNCTION .......... A .......... 19 2.6.1 SQUARED SUM OF RESIDUALS .................... 19 2.6.2 PRIOR INFORMATION REGARDING THE PARAMETERS .. 20 2.6.3 REGULARIZATION .............................. 21 MINIMIZATION OF THE OBJECTIVE FUNCTION ............. 22 2.7.1 GAUSS METHOD OF MINIMIZATION ................ 23 2.7.2 STEEPEST DESCENT METHOD ..................... 25 2.7.3 CONJUGATE GRADIENT METHOD ................... 25 ii 2.8 EXTENSION TO GENERAL FUNCTION ESTIMATION PROBLEM ... 26 CHAPTER 3 CALCULATION OF SENSITIVITY COEFFICIENT AND GRADIENT OF OBJECTIVE FUNCTION ................................... 28 3.1 INTRODUCTION ....................................... 28 3.2 PHYSICAL MEANINGS OF SENSITIVITY COEFFICIENT AND GRADIENT .............................................. 30 3.3 DIRECT METHOD ...................................... 31 3.4 ADJOINT METHOD ..................................... 33 3.4.1 DOGRU AND SEINFELD'S APPROACH ............... 34 3.4.2 HAUG AND ARORA'S APPROACH ................... 38 3.4.3 ALIFANOV'S APPROACH ......................... 40 3.5 FINITE DIFFERENCE METHOD ........................... 44 3.6 CONCLUSION ......................................... 45 CHAPTER 4 IHCP2D: COMPUTER PROGRAM FOR SOLUTION OF GENERAL TWO-DIMENSIONAL INVERSE HEAT CONDUCTION PROBLEMS ........ 47 4.1 INTRODUCTION ....................................... 47 4.2 INVERSE HEAT CONDUCTION PROCEDURE .................. 49 4.3 USING EXISTING DIRECT PROBLEM SOLVER ............... 59 4.4 INPUT/OUTPUT FILES ................................. 60 4.4.1 INPUT VARIABLES OF INIHCP.DAT ............... 61 4.4.2 INPUT VARIABLES OF INTPZ.DAT ................ 65 4.4.3 OUTPUT FILE: OUTIHCP.DAT .................... 65 4.5 DISCUSSIONS ON SOME INPUT VARIABLES OF INIHCP.DAT .. 67 4.5.1 MINIMAL TIME AND MAXIMAL TIME ............... 68 4.5.2 NUMBER OF FUTURE TIME STEPS ................. 69 4.5.3 RATIO OF TIME STEP SIZES .................... 70 iii 4.5.4 MAXIMUM NUMBER OF ITERATIONS ................ 71 4.5.5 REGULARIZATION CONSTANTS .................... 71 4.5.6 LOCATIONS OF SENSORS ........................ 75 4.5.7 WEIGHTING CONSTANTS OF SENSORS .............. 75 4.6 EXAMPLE PROBLEM: TWO-DIMENSIONAL LEADING EDGE GEOMETRY .......................................... 77 4.7 CONCLUSIONS ........................................ 95 CHAPTER 5 SOLUTION OF ONE-DIMENSIONAL INVERSE MELTING PROBLEM 97 5.1 INTRODUCTION ....................................... 97 5.2 PROBLEM DESCRIPTION ................................ 98 5.3 SEQUENTIAL FUNCTION SPECIFICATION METHOD ........... 102 5.4 IMPLEMENTATION OF MELT TO PROGRAM IHCP2D ........... 104 5.5 SENSITIVITY COEFFICIENTS ........................... 105 5.6 RESULTS AND DISCUSSIONS ............................ 110 5.7 SUMMARY AND CONCLUSIONS ............................ 117 CHAPTER 6 CONCLUSIONS AND FUTURE STUDIES .................... 118 6.1 CONCLUSIONS ........................................ 118 6.2 FUTURE STUDIES ..................................... 122 APPENDIX A SUBROUTINE HIERARCHY ............................. 124 APPENDIX B INTERFACE TRACING METHOD FOR THE SOLUTION OF ONE-DIMENSIONAL PHASE CHANGE PROBLEM ................... 128 8.1 INTRODUCTION ....................................... 128 B.2 PROBLEM DESCRIPTION ................................ 130 3.3 INTERFACE TRACING METHOD ........................... 131 B.4 DERIVATION OF FINITE DIFFERENCE EQUATION FOR INTERFACE NODE ........................................ 133 iv B.5 EXAMPLE PROBLEM .................................... REFERENCES LIST OF TABLES Table 4.1 Summary of Input Variables for INIHCP.DAT ......... 66 Table 4.2 Values of Nf, ”0’ and ”l for Stable Algorithm ..... 74 Table 4.3 Material Properties and Some Constants Used in Example Problem ................................... 80 Table 5.1 Summary of the Investigation Results .............. 115 Table 3.1 Finite Difference Equations for all the Nodes ..... 136 vi Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure LIST OF FIGURES 2.1 Geometry of General Two-Dimensional Inverse Problem .......................................... 4.1 Flowchart for the Sequential Estimation of Parameters ....................................... 4.2 Geometry and Nodal Definitions of Example Problem 4.3 Effects of Number of Future Time Steps on Estimated Heat Flux at Node 23 ................... 4.4 Effects of Number of Future Time Steps on Estimated Heat Flux at Node 2 .................... 4.5 Effects of Number of Future Time Steps on Estimated Heat Flux at Node 19 ................... 4.6 Effects of Number of Iterations on Estimated Heat Flux at Node 23 .................................. 4.7 Effects of Number of Iterations on Estimated Heat Flux at Node 2 ................................... 4.8 Effects of Number of Iterations on Estimated Heat Flux at Node l9 .................................. 4.9 Effects on Estimated Heat Fluxes due to Single Error in Sensor 1 at Time - l s .................. 4.10 Effects on Estimated Heat Fluxes due to Single Error in Sensor 2 at Time - 1 s .................. 4.11 Effects on Estimated Heat Fluxes due to Single Error in Sensor 3 at Time - l s .................. 5.1 Geometry of One-Dimensional Inverse Melting vii Figure Figure Figure 5.4 5.2 5.3 Figure A.l Figure Figure 8.1 B.2 Problem .......................................... 101 Sensitivity Coefficient Histories of u ........... 108 Sensitivity Coefficient Histories of T(xs,t) ..... 109 Comparison of Estimated Results with True Values . 114 Subroutine Hierarchy of Program IHCP2D ........... 125 Geometry of One-Dimensional Melting Problem ...... 129 Comparison of Numerical and Analytical Results ... 139 viii 11 31k a1, a2, a3 N) f() NOMENCLATURE 1225211113213 coefficient matrix in equation (4.16) area entry of matrix A slab thickness of one-dimensional problem polynomial coefficient filter coefficients vector of initial guesses of parameters or prior information of parameters right hand side vector in equation (4.16) entry of vector C specific heat marching direction bias expected value function heat generation matrix related to the first order regularization matrix related to the second order regularization heat transfer coefficient identity matrix ix thermal conductivity Lagrange polynomial latent heat number of subdomains number of future measurements used to estimate present parameters degree of polynomial for the jth subdomain number of parameters number of sensors number of time steps outward normal vector on the surface x component of n y component of n iteration number heat flux total length along the jth subdomain local surface coordinate for the jth subdomain the kth location of which the nodal heat flux has been defined as parameter objective function vector stared the calculated temperatures initial temperature calculated temperature of the jth sensor at time ti time final measurement time At m At c measurement time step size calculation time step size weighting coefficient matrix for prior information internal energy volume variance weighting coefficient matrix for measurements weighting constant for the jth sensor sensitivity coefficient of the jth sensor temperature at time t1 with respect to the kth parameter x-coordinate sensor location for one-dimensional problem y-coordinate measurement from the jth sensor taken at time ti thermal diffusivity the kth parameter at time tm vector of estimated parameter values vector of parameter values at the nth iteration small perturbation in the parameter change in the parameter from one time to another variance-covariance matrix of the errors marching step size measurement error corresponding to measurement Y.i emissivity of the surface adjoint variable temperature or interface location xi the zeroth order regularization constant the first order regularization constant interface location density standard deviation or Stefan-Boltzmann constant estimated standard deviation Dirac delta function time mean squared error boundary spatial domain xii CHAPTER 1 INTRODUCTION A conventional heat conduction problem is to determine the interior temperature distribution of a solid with prescribed boundary conditions. An inverse heat conduction problem (IHCP), conversely, is one utilizing transient temperature measurements inside the solid or on the surface to calculate the unknown boundary conditions. The conventional heat conduction problem often is referred to as the direct problem to be distinguished from the inverse problem. The purpose of this dissertation is to present solution methods for general two- dimensional inverse heat conduction problems and one-dimensional inverse melting problems. The following describes the outline of this chapter. Section 1.1 discusses the applications of the inverse heat conduction problem. Section 1.2 presents a literature review for the solution methods of inverse heat conduction problems, and the last section outlines each chapter of this dissertation. 1.1 APPLICATIONS An inverse heat conduction problem arises from the calculation of the surface heat flux on reentering vehicles. Another related problem is the estimation of the impinging heat flux on the leading edge of a hypersonic aircraft. Since the temperatures on the leading edge can be as high as 1500 X, the effects of ionization of the gas molecules and heat radiation become important. These effects coupled with turbulent flow around the leading edge further complicate the problem which is already difficult to solve. An indirect approach is to calculate the 'effective' surface heat flux irrespective of what actually happens. Since the thermocouples will not survive under this harsh environment, they usually are embedded inside the solid. Given these temperature measurements, the inverse heat conduction procedure is used to find the heat flux. Some studies related to this application were published by Mulholland et a1. (1975), Villiams and Curry (1977), Imber (1979), and Blackwell (1981). The second application is the moulding process. Often, it is desired to have the specific movement or propagation rate of the phase change interface during moulding. This is an inverse design problem, and proper boundary conditions are to be found to achieve this specific condition. A similar procedure can also be applied to electronics cooling. Due to highly intense energy generated by electronic components, the removal of this excessive energy is crucial to the performance of components. In other words, boundary conditions are to be found for energy extracting processes. The inverse heat conduction procedure can also be applied to estimate the heat transfer coefficient. Tryanin (1987) and Alifanov et a1. (1987) have studied the heat transfer coefficients of a porous body. In another paper, Osman and Beck (1987) estimated the heat transfer coefficients in the quenching process of a copper ball based on transient temperature measurements inside the body. Since the unknown boundary condition generally is a function of one or more variables, the same procedure can also be employed to solve function estimation problems. Lee et al. (1986) have estimated the spatial distribution of permeability and porosity of an areal petroleum reservoir based on the pressure history measurements at several wells. 1.2 LITERATURE REVIEW Several techniques are available for the solution of IHCPs. In the exact matching method proposed by Stolz (1960) , unknown boundary conditions are solved by setting the measured temperatures equal to the calculated temperatures. This algorithm is unstable as time steps are made small and is very sensitive to temperature measurement errors. The second technique is the function specification method in which a f“Fictional form is assumed for the unknown heat flux, and the best fit function is sought from a family of candidates. Publications "'Ploying this method include Frank (1963), Beck et a1. (1982), Beck, Blackwell, and St. Clair (pp. 119, 1985.1), Flach and Ozisik (1986), and Osman and Beck (1987) . Another widely used technique is the regularization method in which the unknown function is discretized into many components and these components are estimated simultaneously. Many authors have reported their studies on this approach including Alifanov (1973), Alifanov and Artyukhin (1975), Tikhonov and Arsenin (1977), Romanovskii (1984), Beck, Blackwell, and St. Clair (pp. 134, 1985.1), Beck et a1. (1985.2), Alifanov (1986), Markin and Pyatyshkin (1987), and Hensel and Hills (1987). .The last technique to be discussed is the trial function method which includes the function specification method and the regularization method. This method incorporates prior information regarding the functional form of the estimated heat flux. Related works have been reported by Twomey (1963 and 1965) and Beck, Blackwell, and St. Clair (pp. 145, 1985.1). The first part of this dissertation concentrates on the solution of two-dimensional inverse heat conduction problems; publications of this type of problem includes Imber (1974 and 1975), Bass and Ott (1980), Alifanov and Kerov (1981), Lazuchenkov and Shmukin (1981), Kerov (1983), Tang (1985), and Busby and Trujillo (1985). Another important but little studied area of the IHCP is the inverse phase-change problem. The same algorithm solving two- dimensional IHCPs also can be used to solve the inverse phase-change problem if a proper direct problem solver is employed. The methods of solving the direct phase-change problem include the analytic method (pp. 406, Ozisik, 1979), integral method (pp. 416, Ozisik, 1979), moving heat source method (pp. 423, Ozisik, 1979), enthalpy method (Voller and Cross, 1981), heat capacity method (Bonacina et a1., 1973), and interface tracing method (Katz and Rubinsky, 1984, and Goodrich, 1978). Publications on the inverse phase-change problem include Szczurek (1979), and Bobula and Twardowska (1985). 1.3 OUTLINE Four major tasks are to be accomplished in this dissertation. The first one is to examine the function specification method. The function specification method has been employed in several studies to solve one-dimensional inverse heat conduction problems. However, its application to two-dimensional problems is yet to be explored. The second task is to compare different methods of calculating sensitivity coefficients and the gradient of the objective function. The calculations of these quantities are computationally expensive. Therefore, an efficient approach should be used to evaluate these quantities. The next assignment is to develop a general two- dimensional inverse heat conduction code using a combined method of function specification and regularization. Two-dimensional inverse heat conduction problems then can be solved using this general inverse code. The last task is to solve the one-dimensional inverse melting problem using the same inverse code but with a different direct problem solver. The following describes each chapter in this dissertation. In Chapter 2, the function specification method is discussed in detail. The algorithm is derived from these steps: (a) approximating the unknown function with a functional form, (b) estimating the parameters sequentially or simultaneously, (c) defining the objective function, and (d) minimizing the objective function. Different approaches in each step lead to different function specification algorithms. Nevertheless, these algorithms all assume a functional form for the unknown heat flux which is the basic idea of the function W method. Chapter 3 derives and compares three approaches for the calculation of sensitivity coefficient in the inverse heat conduction procedures. These are the direct method (pp. 411, Beck and Arnold, 1977), the adjoint method (Dogru and Seinfeld, 1981, pp. 167 and 335, Haug and Arora, 1979, Sec. 6.5.2, Alifanov, 1986, and Hensel and Hills, 1987), and the finite difference method (pp. 410, Beck and Arnold, 1977, Sec. 6.4, Alifanov, 1986). In Chapter 4, a description is given of a general two-dimensional inverse heat conduction program which is written based on the combined method of function specification and regularization. The spatial distribution of the unknown heat flux is approximated by the piecewise- polynomial function, which can be any combination of different degree polynomials. With the piecewise-polynomial approximation, the discontinuity in heat flux values can also be treated. In general, an inverse problem program is comprised of the inverse heat conduction procedure and a direct problem solver. A useful programming strategy is also described in Chapter 4. The key concept is to use an existing code as the direct problem solver. The general two-dimensional heat conduction code TOPAZ (Shapiro, 1984) is used in this study for this purpose. The resulting inverse problem code, IHCP2D, is capable of solving very general two-dimensional inverse heat conduction problems. Input features of the program are also explained in Chapter 4. With this concept, the function estimation schemes can be used to solve other function estimation problems by using a proper direct problem solver. Chapter 5 investigates a one-dimensional inverse melting problem. The geometry of the problem is a finite slab. An unknown heat flux is applied to the left boundary, and the right boundary is insulated. Both the temperature measurements and the measurements of the interface locations are used to estimate the unknown heat flux. Notice that two different types of measurements have been used to estimate the desired function. Since the interface measurements are used, this procedure can be used to find proper boundary conditions that generate desired interface movement. The interface tracing method (Goodrich, 1978) is employed to solve the direct problem because it traces the interface location at every moment. By replacing TOPAZ with this phase change problem solver, program IHCPZD is modified to solve the inverse melting problem. The effects on the estimated results with respect to several inverse problem parameters, such as the number of future time steps and the measurement time step sizes, are also studied. The last chapter, Chapter 6, concludes this dissertation. Some suggestions for future studies also are presented. These future studies include: the development of a computer program for the solution of three-dimensional inverse heat conduction problems; the optimal design study to find optimal values for locations of sensors, number of sensors, etc.; and employing numerical analysis to analyze the errors, stability, convergence, and consistency of the inverse problem algorithms. CHAPTER 2 GENERAL DISCUSSION OF THE FUNCTION SPECIFICATION METHOD 2.1 INTRODUCTION Four techniques for the solution of inverse heat conduction problems have been discussed by Beck et al. (Chapter 4, 1985.1). These procedures are the single future time step method, the function specification method, the regularization method, and the trial function method. In this chapter, the function specification method is considered for the solution of the general two-dimensional inverse heat conduction problem. The known boundary conditions of this problem include the first kind (prescribed surface temperatures), the second kind (prescribed surface heat flux), and the third kind (convective boundary condition). An unknown heat flux is also applied to part of the boundary. The objective is to estimate this unknown surface heat flux based on some temperature measurements taken from the solid. Although the two-dimensional geometry has been used, the procedure can be applied to three-dimensional problems as well. Many studies have been published on the function specification method including Frank (1963), Beck et a1. (1982), Flach and Ozisik (1986), and Osman and Beck (1987). In this chapter, the derivation of this method is discussed in great detail, and it is found that different variations in each step of the derivation lead to different function specification methods. Later, in Chapter 4, a particular function specification method is introduced and employed to develop a computer program for the solution of general two-dimensional IHCPs. The following describes the outline of each section of the chapter. Section 2.2 describes the problem and presents its mathematical formulation. Section 2.3 outlines the function specification method. Section 2.4 discusses the functional form used to approximate the unknown heat flux. The parameters associated with the approximating function can be estimated sequentially or simultaneously. A comparison of these two approaches of estimation is presented in Section 2.5. Section 2.6 discusses the definition of the objective function. Section 2.7 presents three methods for minimizing the objective function, and the last section summarizes these discussions. 10 Figure 2.1 Geometry of General Two-Dimensional Inverse Problem 11 2.2 PROBLEM DESCRIPTION In this section, the mathematical description of a general two- dimensional inverse heat conduction problem is given. The geometry considered is shown in Figure 2.1. The spatial domain is denoted by 0. The symbols P1, F2, and P3 represent boundaries with known boundary condition of the first kind (prescribed surface temperatures), second kind (prescribed surface heat flux), and third kind (convective boundary condition), respectively. Furthermore, any known radiative boundary condition can be treated as a convective boundary condition by defining a temperature-dependent radiative heat transfer coefficient. An unknown heat flux has been applied to the boundary F4. The initial temperature distribution of the solid is f4(x,y). The objective of this problem is to estimate the unknown surface heat flux on boundary 1‘4 given surface and/or subsurface transient temperature measurements. The mathematical model of this problem is d_, d1 1. ET ax (Rx ax) + 8y (ky 6y ) + g(X.y.t) - pc %% (x,y) in 0, 0 s t s tf (2.1) T(x,y,t) - f1(x,y,t) (x,y) on F1, 0 s t S tf (2.2) - 31 - 31 - kx 3x nx ky 8y ny f2(x,y,t) (x,y) on F2, 0 s t s tf (2.3) - 31 - 111' - - kx 3x nx ky 8y ny h[T(x,y,t) f3(t)] (x,y) on F3, 0 s t s tf (2.4) T(x,y,0) - f4(x,y) (x,y) in 0 (2.5) -kx g—E nx «y 3% ny - q(x,y,t) - 7 (x,y) on r4, 0 s c 5 cf (2.6) 11 2.2 PEBBLE! DESCRIPTION In this section, the mathematical description of a general two- dimensional inverse heat conduction problem is given. The geometry considered is shown in Figure 2.1. The spatial domain is denoted by O. The symbols P1, P2, and P3 represent boundaries with known boundary condition of the first kind (prescribed surface temperatures), second kind (prescribed surface heat flux), and third kind (convective boundary condition), respectively. Furthermore, any known radiative boundary condition can be treated as a convective boundary condition by defining a temperature-dependent radiative heat transfer coefficient. An unknown heat flux has been applied to the boundary P The initial 4. temperature distribution of the solid is f4(x,y). The objective of this problem is to estimate the unknown surface heat flux on boundary F4 given surface and/or subsurface transient temperature measurements. The mathematical model of this problem is El a"(k a1) + a‘ (k ix) + 3(X.y.t) - pc at 3x x ax 8y y 6y (x,y) in 0, 0 S t s tf (2.1) T(x,y,t) - f1(x,y,t) (x,y) on F1, 0 s t S tf (2.2) -fiI-3-T- kx ax nx ky 3y ny f2(x,y,t) (x,y) on F2, 0 s t s tf (2.3) _ i1 _ flI _ - kx ax nx ky 3y ny h[T(x,y,t) f3(t)] (x,y) on P3, 0 s t s tf (2.4) T(X.y.0) - f4(X.y) (X.y) in 0 (2.5) -fl-§1- - kx ax nx ky 6y ny q(x,y,t) ? (x,y) on F4, 0 5 t s tf (2.6) 12 where q(x,y,t) is the unknown heat flux to be estimated. The transient temperature measurements are in - T(xj’yj’ti) + £11 1 s j 5 NS, 1 s i s N (2.7) t where N8 is the number of sensors, NC is the number of time steps, in is the temperature measurement, and 611 is the measurement error. The point (xj time when the measurement is taken. ,yj) is the location of the jth sensor, and t represents the i In the above equations, T represents the 'true' temperature. The quantity g(x,y,t) represents the heat generation, which is a known function of position and time. The symbol tf represents the final time. The thermal properties are thermal conductivity, k, specific heat, c, and heat transfer coefficient, h. These properties, k, c, and h, in general, are known functions of temperature; h can also be a known function of time. The subscripts 'x' and 'y' of k indicate that the material can be orthotropic. Density is denoted by p. The f functions, f f 1, 2, 3, and f4 are known. The symbols nx and ny are the x- and y-components of the outward normal unit vector u. With this definition of n, the heat flux legging a surface is positive. For the clarity of notation, a surface coordinate r is introduced, which traces along the surface F With this surface coordinate, the 4. unknown surface heat flux can be expressed as q(X.y.t) - q(r,t) (x,y) on Pa, O S r S R (2.8) 13 where R represents the total length along the boundary FA' 2.3 OUTLINE OF THE FUNCTION SPECIFICATION METHOD The function specification method is used to estimate the unknown surface heat flux q(r,t) on F4. The following gives an outline of the procedure. More extensive discussions will be presented in the subsequent sections. Step 1: Assume a functional form for the unknown heat flux. The basic concept of the function specification method is to specify the unknown heat flux with an assumed functional form. For example, the spatial variation of the unknown heat flux can be approximated by a second degree polynomial as q - [3fi1(t) - 452 + 33(c>1 fi 2 + 2[pl(c) - 252(c) + 53(c)] :2 o s r s R, o s c s t (2.9) where 51(t) - q(0,t), 52(c) - q(R/2,t), and fi3(t) - q(R,t). Note that this is not the only way the parameters can be defined. Other options such as the polynomial coefficients and derivative values can also be defined as the parameters. Next, the temporal variations of these parameters are specified with another functional form. For the example in equation (2.9), piecewise-constant functions are used and 14 flk(t)-,Bkm lsks3,15msNt,tm_11+(b — a)": u (b - p) 2 T +§ (01(Nifl) nip (2.19) 1-0 where the scalar ”i is the ith regularization constant. The zeroth order regularization has Nb - I, the identity matrix. For a case with Np - S, the first order regularization has the N1 matrix of 22 -1 1 O 0 0 0 -l l 0 0 H1 - 0 0 -1 1 0 (2.20a) O 0 0 -l l 0 0 0 0 0 and the second order regularization has 1 -2 1 0 0 0 l -2 1 0 H2 - 0 0 1 -2 1 (2.20b) 0 0 0 0 0 0 0 0 0 0 In addition to the above options, prior information can be incorporated in the regularization terms rather than a term by itself as shown in equation (2.18). This approach, the trial function method, has been reported by Twomey (1963 and 1965). 2.7 MINIMIZATION OF THE OBJECTIVE FUNCTION Extensive studies have been performed in finding the minimum of an objective function including Haug and Arora (1979) and Alifanov (1986). Results from these studies can be readily employed to minimize the objective function. In general, these methods seek the minimum of the objective function by selecting the proper direction and step size for a piecewise path in the parameter space, which, hopefully, will eventually lead to the minimum. The general procedure is described in the following. Given the initial guess of the parameters 5(0), the minimum is calculated using the iterative scheme p(n+1) _ pm) + 1(“+1)d(“+1) (2.21) 23 where the superscript n stands for the nth iteration. The vector d] v [Y - rm 1+ “(a )5 ]} + u b (2.27) where the superscript n represents the iteration number. This is the result of the Gauss method. Equation (2.27) represents a set of simultaneous equations with p being the unknown vector. Notice that the temperature T and the sensitivity matrix g% are calculated at 3(n). Techniques for solving a set of linear equations can then be employed to find the solution. Among these techniques are the Regula- Falsi method (p.76, Conte and De Boor, 1979), Newton method (p.79, 25 Conte and De Boor, 1980), and fixed point iteration method (p.88, Conte and De Boor, 1980). The solution will give both the marching direction and the step size. 2.7.2 STEEPEST DESCENT.METNOD In the steepest descent method (Section 6.4, Alifanov, 1986 and p.45, Haug and Arora, 1979), the marching direction is chosen to be opposite to the gradient of the objective function and the step size is selected to minimize the objective function along this direction; that is, (n+1) 15. (n) (n+1) _ min (n) _ Qfi (n) ‘ ‘1 .1 S[fl 7 dpw )] 0-”) where d is the marching direction and 7 is the step size. The superscript n stands for iteration number. Equation (2.28) states that the marching direction is in the opposite direction of the gradient, and equation (2.28) implies that the step size, 1, is chosen to be the value that minimizes 8 along the line passing 8(a) in the direction of d. Equation (2.29) represents a one-dimensional minimization problem which searches for the 1 value minimizing the objective function along AS the direction dfl' 2.7.3 CONJUGATE GRADIENT METHOD 26 The marching direction in the conjugate gradient method (Section 6.4, Alifanov, 1986 and p.91, Haug and Arora, 1979) uses a linear combination of the gradient of the objective function and the direction of the previous iteration, in other words, d(n+l) _ w(mndm) _ £5;me (2.30) where "(n+1) - I g%(p:n:l))2|2 |§%w“| (2.31) The scalar w is a weighting constant. The step size, 1, is still calculated using equation (2.29). 2.8 EXTENSION TO GENERAL FUNCTION ESTIMATION PROBLEM The function specification method discussed in this chapter has the potential of solving the general funggign_gggim§§;gn problem. The problem formulation, the unknown function, and the measured quantities might be different, but the algorithm is basically the same. For example, the unknown function might be the surface temperatures, or the measurements might be the velocities for a momentum-energy coupled problem or the interface history of a phase-change problem. The exploitation of the application of this method to many other function estimation problems is a challenging and fruitful prospect. 2.9 SUMMARY 27 A summary of the function specification procedure is the following. 1. Assume a functional form for the unknown heat flux. 1.1 Decide whether global or subdomain specification is to be used. 1.2 Decide what type of functions are to be used. 2. Employ whole domain or sequential estimation. 3. Define the objective function. 3.1 Define the deviation measure. 3.2 If there is any prior information regarding the parameters, include it in the objective function. 3.3 Decide which regularization procedure, if any, is to be used. 4. Minimize the objective function. 5. Implement the procedure in a computer program. Variations in each step lead to different function specification methods. CHAPTER 3 CALCULATION OF SENSITIVITY COEFFICIENT AND GRADIENT OF THE OBJECTIVE1FUNCTION 3.1 Imopucnon In Section 2.7, three methods were presented to minimize the objective function. The first one is the Gauss method which requires the calculation of the sensitivity coefficient, aTi/apj. The other two are the steepest descent and the conjugate gradient methods which require the calculation of the gradient of the objective function, BS/OBJ. Information provided by the sensitivity coefficient and the gradient carries significant physical insight which will be discussed in Section 3.2. Several techniques have been developed to calculate the sensitivity coefficient and the gradient. The first technique is the direct method which calculates the sensitivity coefficient by solving the sen:1;121§1_pxghlgm. The gradient then is found from a relation 28 29 between sensitivity coefficient and the gradient. This method will be described in Section 3.3. The second technique is the adjoint method which solves for the sensitivity coefficient or the gradient through an adjoint variable, A, which is the solution of the ggjgin;_pxg§lgm. Alifanov (1986) has employed this method to calculate the gradient for the estimation of parameter function in the function space. Haug and Arora (1979) also have applied this method to calculate the gradient to solve optimization problems in mechanical systems (p.167 and p.335, Haug and Arora, 1979). Dogru and Seinfeld (1981) used this method to calculate the sensitivity coefficient rather than the gradient. Hensel (1986) also has applied this method to calculate the sensitivity coefficient. Some of these techniques will be described in Section 3.4. The last technique to be discussed is the finite difference method which approximates the partial derivative by the finite difference equation. This method will be presented in Section 3.5. Throughout this chapter, a one-dimensional transient linear problem is used as an example. The same procedures can be applied to more complicated problems as well. The mathematical formulation of the problem is 231: :11 a a 2 - at 0 < x < a, 0 s t s tf (3.1) x T(x,0) - T0(x) 0 s x S a (3.2) 4: 330,12) - q(t;p) 0 s t 5 Cf (3.3) 30 -k %§(a,t) - qa(t) o s c 5 cf (3.4) The heat flux applied to the surface x - 0 is unknown and depends on the parameter vector p'- [61], l S i s Np. The heat flux applied to the surface x - a is known and does not depend on the parameters. To be distinguished from the sensitivity problem and the adjoint problem to be discussed, equations (3.l)-(3.4) is called the temperature Hrshlem. The objective function is s - I [Y(t) - T(xs,t)]2dt (3.5) t where x8 is the sensor location. The transient temperature measurements at this sensor location are denoted by Y(t). This definition can be modified to include multiple sensors. For multiple sensors, the subscript s is from 1 to N3, the number of sensors. 3.2 PHYSICAL.MEANINCS OF SENSITIVITY COEFFICIENT AND GRADIENT The sensitivity coefficient was defined in equation (2.24) as ari/apj, which represents the changes in T1 with respect to the changes in the parameter, fij‘ Information provided by the sensitivity coefficient is very useful. Generally, it is preferable to have large, uncorrelated sensitivity coefficients. In other words, sensors should be placed at locations that are most sensitive to the parameters. The 31 sensitivity coefficient also can be used to determine the measurement time step size and the variance of estimated values (p.248, Beck and Arnold, 1977). As pointed out by Beck and Arnold (pp. 349, 1977), if the sensitivity coefficients are linearly dependent, the objective function for the ordinary least square estimation or the maximum likelihood estimation has no unique minimum point. In this case, it is impossible to achieve convergence to unique parameter values. The gradient of the objective function is defined as BS/BBj. This quantity represents the change in the objective function, S, with respect to the parameter, flj' In other words, it indicates that effect of the parameter on the objective function. 3.3 DIRECT METHOD The procedure of the direct method is described in the following. The first step is to differentiate the temperature problem with respect to each parameter, fij. This results in the sensitivity problem, and its solution is the sensitivity coefficient, aT/fij. To calculate the gradient, a relation between the sensitivity coefficient and the gradient is employed. Differentiated with respect to £1, 1 s j s Np, equations (3.1)- (3.4) become 32 02:332-ngng 0 t and the solution A(x,r;tf) is available, then the solution A(x,7;t) can be obtained by coordinate transformation as A(x,r;t) - A(x,r-t +t;tf) (3.24) f With this relation, then actually only one adjoint problem needs to be solved which is A(x,r;tf). Substituting equation (3.24) into equation (3.23), it becomes 37 t :21 - 3510‘s.!» Jo , .2831. A(0,1 tf+t,tf) k a5.10.5) d1 (3.25) This is the equation used to calculate the sensitivity coefficient from the adjoint variable. The adjoint method presented in this section can be summarized in two steps. First, the adjoint variable, A(x,1), is solved from equations (3.18)-(3.21) with t - t Next, the sensitivity f. coefficients are calculated using equations (3.25). Note that if there are Ns sensors each at different locations x8, then there will be Ns adjoint problems, one for each sensor. For an inverse problem with one sensor and Np parameters, there is only one adjoint problem to be solved compared to Np sensitivity problems to be solved using the direct method. To see the physical meaning of the adjoint variable, a coordinate transform, 1' - t - 1, is introduced into equations (3.18)-(3.21). They then become au-fl, 0 -k §§I*511(a.t) - qa Ns’ the adjoint method is preferable. If Np s Ns' the direct method or the finite difference method is preferable. For nonlinear temperature problems and Np 5 N8, the direct method is preferable because it only involves the solutions of Ns linear sensitivity problems and one ‘. - 46 nonlinear temperature problem compared to Ns+1 nonlinear temperature problems if finite difference method were to be used. However, extra programming effort is required for the direct method to solve the sensitivity problem. Three approaches of the adjoint method have been presented. Although similar procedures are used to derive these techniques, the adjoint problems are defined differently. This indicates that there is no unique definition of the adjoint problem, and any definition is justified as long as it can simplify the integral equations (3.17) and (3.54). CHAPTER 4 IHCP2D: COMPUTER PROGRAM FOR.SOLUTION OF GENERAL.THO-DIMENSIONAL.INVERSE HEAT CONDUCTION PROBLEMS 4.1 INTRODUCTION Two popular procedures for the solution of inverse heat conduction problems are the function specification method and the regularization method. The function specification method emphasizes the approximation of the unknown heat flux with a functional form (Section 2.4), and the regularization method involves the implementation of regularization terms to the objective function. Since these two procedures emphasize different aspects of the technique, a function specification method can also be a regularization method if the regularization terms are included in the objective function (Section 2.6) and vice versa. In this chapter, a combined method of function specification and regularization is proposed for the solution of general two-dimensional inverse heat conduction problems. As part of this combined method, a 47 48 piecewise-polynomial function, which can be the piecewise combination of different degree polynomials, is used to approximate the spatial distribution of the unknown heat flux. This combined method has been employed to develop a computer program, IHCPZD, to solve general two-dimensional inverse heat conduction problems. In this program, a special programming concept is adopted. The strategy is to modify an existing and widely-accepted direct problem solver into a subroutine and to incorporate it in the inverse heat conduction algorithm. The general two-dimensional heat conduction code TOPAZ (Shapiro, 1984) is used as the direct problem solver in program IHCPZD. With this powerful direct problem solver, the inverse code IHCPZD is capable of solving a wide variety of practical inverse heat conduction problems. This procedure gives greater power, generality, and reliability with less programming effort than if a two-dimensional direct problem code were developed specifically for the inverse heat conduction problems. This idea of using an existing code has been employed by Hills (1987) to develop a parameter estimation computer program, ESTIM. Nevertheless, it is used for the first time in the area of inverse heat conduction problems in the development of IHCPZD. The contents of this chapter is briefly outlined below. The general two-dimensional inverse heat conduction problem and its mathematical formulation have already been described in Section 2.2. Section 4.2 presents the proposed combined method to solve this problem. The concept of using an existing direct problem code is described in Section 4.3. Section 4.4 explains the input and the output files accessed by code IHCPZD. Some of the input variables for 49 the inverse problem input file are discussed in Section 4.5. Program IHCPZD has been successfully used to solve several inverse problems; two of these examples are included in Section 4.6. The last section gives some conclusions. 4.2 INVERSE HEAT CONDUCTION PROCEDURE This section discusses the proposed combined method of function specification and regularization. The first step is to approximate the spatial distribution of the unknown heat flux with a piecewise-polynomial function. The domain PA, with unknown heat flux, is divided into several subdomains denoted by P P Fa N , where Nd is the total number of subdomains. For 4,1 d 4,2’ "’ the clarity of notation and programming, a local surface coordinate r J is defined for each subdomain F41. This coordinate traces along the boundary raj. The unknown surface heat flux then can be expressed in terms of r as J q(x,y,t) - q(rj,t) (x,y) on F 0 s r s R (4.1) The quantity R.J represents the total length along the boundary raj. The heat flux over each subdomain is approximated by a polynomial which is 50 +1 N q(rj,t) -§. 211151013"1 0 s rJ s Rj, 0 s t .<. tf (4.2) -1 where ajk's are the polynomial coefficients and N1 is the degree of this polynomial. The second step is to define the parameters to be estimated. The parameters associated with each polynomial are defined as the heat flux values at several locations, that is, fijk(t) - q(rjk,t) l S j s Nd, l S k S Nj+l (4.3) Note that an NJth degree polynomial requires the evaluation of N +1 .1 parameters. The polynomial given in equation (4.2) can be expressed in terms of parameters, fljk(t), using a Lagrange formula (p. 38, Conte and de Boor, 1980) as +1 N q(rj.t) -j fljk(t)lk(rj) 0 k-l IA '1 j 5 RJ' 0 S t 5 cf (4.4) where £k(rj) is an thh degree Lagrange polynomial which is defined by n+1 j r - r 2(r)-n J—J-L lsjsN,1skSN+l (4.5) kj 1-1rjk'rji d j ivk 51 This idea of using a general function to approximate the unknown heat flux also has been reported by Bass et a1. (1980). If the function values at the junction of two subdomains are continuous, the number of parameters can be reduced. For either case, the double- subscript parameters, fljk(t), can be re-numbered to become single- subscript parameter, flk(t), l s k s Np, where Np is the total number of parameters for domain F4. Thus, by re-numbering the parameters the derivation that follows will be similar to the sequential regularization method (pp. 141, Beck et a1., 1985). Next, the time dependent arameters, 6 (t), are discretized in P k time by pk(c) - 5km ‘ 1 s k s Np, cm_1 < c s cm (4.6) The objective then becomes the estimation of the parameter constants, fikm's' These parameters are estimated simultaneously in the 'k' index and sequentially in the 'm' index. Assuming that the parameters Bkl’ fik2’ .., pk,m-l’ l s k s Np’ have been estimated, the task now is to estimate the parameters at the mth time step, pkm’ 1 S k S Np' To improve the stability of the algorithm, the combined method uses several future measurements for the estimation of these parameters. The objective function for the estimation of the parameters, flkm's’ is 52 N8 Nf Np 2 2 S '7 "j E (Yj,m+i-1 ' Tj,m+i-1) + ”0} 5km 1-1 1-1 k—l 2 (4.7) N -1 p I ”1 E (flk+l,m ' flkm) k—l where w is the weighting constant for the jth sensor, “0 and wl are .1 the zeroth and the first order spatial regularization constants, respectively; and the integer Nf is the number of future temperature measurements. By adding the zeroth and the first order regularizations to the objective function, this technique combines regularization with function specification. If a large number of zeroth degree polynomials are used and only few sensors are available, it becomes the regularization method. On the other hand, if both regularization constants are set equal to zero, it becomes the function specification method. A paper published by Beck and Murio (1986) presents similar concepts. In their study, regularization was applied to time variation of heat flux while in this procedure the regularization is applied to the gpggigl variation of heat flux. The first term on the right hand side of equation (4.7) has the effect of reducing the differences between measured and calculated temperatures in the least squares sense. The second term is for the zeroth order regularization which has the effect of suppressing oscillation in the estimated values. The third term is the first order regularization which reduces the differences in subsequent parameters. 53 To calculate the sensor temperatures, the parameters are assumed to be temporarily constant in these Nf future time steps: 1 S k S Np (4.8) pkm ' pk,m+1 ' °°° ‘ pk,m+N -1 f (Recall that these fi's represent parameters associated with the spatial variation of the heat flux at time tm') The last step is to minimize the objective function using the Gauss method of minimization (pp. 340, Beck and Arnold, 1977). Differentiating equation (4.7) with respect to each parameter and setting the results equal to zero yields Ns Nf as 2 an,“ ' ' 2 2": 2 “Lam-1 ' Tim-1’ ‘11! 3-1 1-1 Np-l as as +26} ()8 -B)(-1‘-+-L~“--m) 1 k-l k+1,m km as,“ 382m n9 55 + 260 E pkm 55f: - 0 1 s 2 s Np (4.9) k-l where the sensitivity coefficient, X111, is defined as 8T + -1 X - 4‘14— (4.10) 312 352m 54 It is calculated using the following finite difference approximation in program IHCPZD 8T T (.3 4' A5) ' T (l3 ) 5311111111 ,3 4m A6 1“ (4.11) in Note that only the parameter associated with 1 has been perturbed. Suppose that the initial guesses of parameters, b - [bkm]’ l s k s N P are 'close' to the solution of equation (4.9), p - [6km], 1 s k s Np. The temperature, Tj,m+i-1’ and the sensitivity coefficient, ink’ can be approximated by N A p A TJ ,M1-1(p) s: TJ,ID+I-1(b) '1’} xjik(b) (Bk!!! " bkm) (4-12) k-l and Substituting equations (4.12) and (4.13) into equation (4.9) results in Ns Nf Np ' E “j E [éj,m+i-l ' Tj,m+i-l(b) ’ E xjik(b) (5km ' bkm)] inz(b) 1-1 1-1 k-l - wlfll-l,m + (wo + 2w1)fl£m - wlfl£+l,m - 0 l S 2 S Np (4.14) 55 A Unfortunately, b is not always close to p, and iterations are thus A required to solve for n. (For problems which are linear in the parameters, iteration is not needed.) By replacing b with fl(n) - (n+l)_ (n+l)] [6(“)], 1 s k s Np, and 3 with p -[p l S k s Np, equation (4.14) becomes Np Ns (n) (n) (n+1)_ (n+1) (n+1) 2 [28" "13 xjik x111] 5w. 1152 + (“’0 + 2H1)52 k—l 3-1 1-1 _ N p (n+1) (n) (n) (n) ' “151+1,m' Es" 35f [ Yj,m+1-1 ‘ Tj,m+i-l I §( xjik 5m ] x312 1-11-1 k-l l S l S Np (4.15) where 'n' is the iteration index. In this equation, the temperatures (n) (n). Tj’m+1_1's and the sensitivity coefficients ink s are evaluated at ”(n). This equation can be expressed in matrix form as (A.+ 601 + 6131791) p‘“+1) - c (4.16) where ”0 and ”l are the zeroth and the first order regularization constants, respectively. The coefficient matrix A.is A.- [A1k] where A),k - [ Esw wj Sf x§?; xggi ] (4.17) 3-1 1-1 56 and the right hand side vector C is c - [C2] where N N N s f P (n) (n) (n) (n) c! 'E "1 E [Emu-1 ' Tin-+14 *E x111: ”km ] x112 (“'18) 3-1 1-1 k—l where the coefficient matrix A.generally is a full matrix. The identity matrix I and the matrix “1T“1 are 1 o . . o o 1 o . . 1 - . . . . . (4.19) o 1 o o . o 1 ’ 1 -1 o o ‘ -1 2 -1 o . o I o -1 2 -1 o . . a1 31 - . . . . . . . (4.20) . . . . . . o o o -1 2 -1 . o o -1 1 . Equation (4.16) represents a system of Np simultaneous linear equations with p as unknowns. This system of equations can be solved using Gauss-Jordan elimination (p.272, Carnahan et a1., 1969). For a problem that is nonlinear in terms of heat flux components, iteration is needed to find the converged p(“+1). 57 l fll GIVENfig),lsksz l (n) l CALCULATE 'rj’ 1-109] ) I i (n) 7 iii FOR k - l...Np, PERTURB a] CALCULATET l (5(“)+AB) l j,m+i-l km 1 T (5 +63) ' T (5 ) ._.{ CALCULATE xjik z 4.5M Afl 1MJ Iconsmucr A SYSTEM or EQUATIONSI l r SOLVE FOR ,3”, 151:sz I CONVERGE ? NO YES Figure 4.1 flowchart for the Sequential Estimation of Parameters 58 In the above derivations, if both regularization constants, mo and ”1’ are set equal to zero, the procedure is the sequential function specification method. 0n the other hand, if the number of parameters is greater than the number of sensors, at least one of the regularization constants can not vanish, and the procedure becomes the regularization method. Thus, this method conveniently combines regularization with function specification into one formula. Another powerful feature of this method is the use of the sequential estimation procedure which is computationally efficient compared to whole domain (or simultaneous) estimation over all times. Figure 4.1 shows the flowchart for the sequential estimation of the parameters, fikm’ l s k s Np' The following steps describe the procedure. (1) Assign initial guesses to the parameters. The same values are also assigned to the parameters of the next Nf-l time steps. (2) Calculate the temperatures Tj m+i~l with unperturbed parameter values. (3) For k - l ...Np, perturb the parameter, flég), by adding Afi, and calculate the temperatures Tj,m+i-l' (4) Calculate the sensitivity coefficients ink using the finite difference approximation in equation (4.11). (5) Construct the matrix equation according to equation (4.15). (6) Solve for p‘“+1). 59 (7) Check convergence of the parameters. If at least one parameter has (n+1) as initial guesses, and go to step 1. If all not converged, use p parameters have converged, increase the time index m by l, and repeat the procedure to estimate the parameters of the next time step. In general, iteration is not needed to estimate the heat flux components (even for problems with temperature dependent properties), unless radiation boundary conditions are involved. 4.3 ‘USINC EXISTING DIRECT PROBLEM SOLVER The solution of equation (4.15) requires calculated values of the (n) (n) j,m+i-l x temperature T jik and the sensitivity coefficient at a given set of parameters. A direct problem solver subroutine is required for the evaluations of the above quantities. Great effort is required to develop a direct problem solver that can handle general two-dimensional heat conduction problems. Rather than developing a new direct two- dimensional solver, a much more efficient approach is to modify an existing direct problem solver and to use it as a subroutine for the inverse heat conduction algorithm. A required feature of the direct problem solver when used as a subroutine for the inverse heat conduction problems is that the boundary conditions of prescribed temperatures, heat flux, etc., can be arbitrary functions of both space and time. In this work, the general two-dimensional finite element heat conduction code TOPAZ (VAX version, January 1986) is employed as the direct problem solver. The computer program TOPAZ was developed by A.B. Shapiro (1984) of Lawrence 60 Livermore National Laboratory. TOPAZ can solve the problems with (a) arbitrary two-dimensional geometries, (b) variable temperature boundary conditions, (c) variable flux boundary conditions, (d) variable convective boundary conditions, (e) radiative boundary conditions, (f) temperature dependent properties, (g) composite materials, (h) orthotropic materials, (1) internal contact resistances, and (j) chemical reactions. Advantages of using existing codes are that they can be very powerful, general and reliable, and have been widely accepted and used. Furthermore, the user's manual for a program, such as TOPAZ, has been refined through interaction with many users and also many engineers are familiar with the use of TOPAZ. As a result, the inverse problem program is capable of solving problems with the same complexity as for the direct code. Appendix A presents the subroutine hierarchy of program IHCPZD. In order to use an existing program as a subroutine, some changes are needed. However, to insure correctness of the program, minimal changes are required, and the subroutine version of TOPAZ is called SUBTPZ. A similar procedure can be used to utilize other direct problem solvers to perform function or parameter estimation in diverse areas such as combustion and fluid mechanics. 4.4 INPUT/OUTPUT FILES Two input files are required to execute IHCP2D, The first input file is called INIHCP.DAT which contains the input data for the inverse heat conduction algorithm. All data but the title are read with £§g§_ formats. The second one is INTPZ.DAT which contains the information needed by the direct problem solver SUBTPZ; the input variables and 61 formats of INTPZ.DAT are exactly the same as those defined in pp. 39-59 of TOPAZ user's manual (Shapiro, 1984). The following presents the outline of this section. The input variables of the input file INIHCP.DAT are explained in Section 4.4.1. Since the input variables of the input file INTPZ.DAT are exactly the same as that of program TOPAZ, the TOPAZ user's manual should be used to create the input file INTPZ.DAT. Next, the output file of program IHCPZD is discussed. The last subsection describes the compiling and execution procedure of program IHCPZD. 4.4.1 INPUT VARIABLES OP IRIBCP.DKT Each input variable in the input file INIHCP.DAT is presented in this subsection. There are seven basic blocks of data labeled 0 to 6. The possible ranges of some input variables are enclosed by parentheses, and a summary of these input variables is given in Table h.1. filQ§k_Q: Case Title of the Inverse Problem TITLE TITLE: Reader to appear in output file which can be any letter or number and can include up to 72 spaces on one line. £12£k_1: Control Data THIN, IHAX, NP, NCALn HAXITR, IPNITR THIN: Minimal Time (2 O) TMAX: Maximal Time (5 cf, the final measurement time) 62 NF: Number of future measurements used in sequential estimation of heat flux (Section 4.2). Typical values for NF are between 2 and 5. (1 5 NF 5 10) At NCAL: NCAL - XE“ where Atm is the measurement time step size, c and Atc is the calculation time step size used in the direct problem solver. The purpose of NCAL is to give the Atc value. (NCAL z 1) MAXITR: Maximal number of iterations allowed for the convergence of parameters. For linear problem, set MAXITR - O. For nonlinear problems, suggested value of MAXITR is 2. (MAXITR 2 O) IPNITR: Flag for printing iteration progress of parameters - 1: Print iteration process to the output file. - 0: Do not print. filggk_2: Regularization Constants OHEGAO, OMEGAI OMEGAO: The zeroth order regularization constant, mo (2 0) OHEGAl: The first order regularization constant, ml (2 0) mm: Piecewise-Polynomial Approximation for Unknown Heat Flux nsun mesa), ma) unsun(1,1), unsung), unsun(1,m(1)) NDBETA(1,1), NDBETA(1,2), unam(1,msc(1)+1) Inscmsun), mums“) RDSUB(NSUB,1), NDSUB(NSUB,2), ..., NDSUB(NSUB,NND(NSUB)) 63 NDBETA(NSUB,1), NDBETA(NSUB,2), ..., NDBETA(NSUB,IDEG(NSUB)+1) NSUB: Number of subdomains (l s NSUB s 50) IDEG(I): Degree of polynomial for the Ith subdomain (0 s IDEG s 4) NND(I): Number of nodes on the Ith subdomain (1 s NND 5 100) NDSUB(I,J): Nodal number of the Jth node on the Ith subdomain NDBETA(I,J): Nodal number of the Jth parameter on the Ith subdomain Note: The nodal numbers should be specified in counterclockwise order according to the nodal definitions of the geometry. Blggk_fi: Filtering Information IP23, P281, P232, PRE3 IPOST, POSTl, POSTZ, POST3 IPRE: Flag for pre-filtering process - 1: Perform filtering process on the measurement data. - 0: No filtering process PREI, PRE2, PRE3: Pro-filter coefficients IPOST: Flag for post-filtering process - 1: Perform filtering process on the estimated results - O: No filtering process POSTl, POSTZ, POSTS: Post-filter coefficients The pre-filtering process is defined as Y1, - PREl-YLi-l + PREZ-in + Panza-Y1.“1 (4.21) where ? stands for the filtered quantity. Similarly, the post- filtering process is defined as 64 - POSTl-fik m- + POSTZ-flkm + P081306k m+l (4.22) 3km 1 where 3 stands for the filtered parameter value. In both cases, the sum of filter coefficients should be equal to one. nlggkgi: Sensor Information NS XS(1), YS(1). “(1) 28(88), YS(NS), U(NS) NS: Number of sensors (1 5 NS 5 20) XS(J): X-coordinate of the Jth sensor YS(J): Y-coordinate of the Jth sensor The sensor can be placed anywhere in the body, but two sensors can not be at the same location. W(J): Weighting constant for the Jth sensor. Suggested value is the reciprocal of the variance of the measurement error. If the variances of measurement errors are not known, just enter 1 for all sensors to bypass this feature. filggk_§: Temperature Measurements TIME(1), Y(1,1), ...., Y(NS,1) TIHE(NT),‘Y(1,NT), ...., Y(NS,NT) TIME(I): Time t1 Y(J,I): Temperature measurement taken from the Jth sensor at time ti 65 4.4.2 INPUT VARIABLES 0P INTPZ.DAT The input file INTPZ.DAT contains the information required by subroutine SUBTPZ. Since SUBTPZ is a subroutine version of TOPAZ, the input variables of INTPZ.DAT are exactly the same as those for TOPAZ. A detailed description of the input variables for TOPAZ is given in the TOPAZ user's manual (Shapiro, 1984). Although TOPAZ can solve problems with either uniform or variable time step sizes, program IHCPZD is restricted to uniform time step sizes over the time domain. 4.4.3 OUTPUT FILE: OUTIHCP.DAT The output file is divided into three major sections. The first one records the input variables from the input file INTPZ.DAT. The second part lists the input variables from the other input file INIMCP.DAT. The last part contains the estimated results including estimated parameters, estimated surface nodal heat fluxes, and comparison of the measured and calculated temperatures. The temperature residuals and the root-mean-squares of the residuals are also calculated and printed at the end. 66 Table 4.1 Summary of Input Variables for INIHCP.DAT filggk_9: (Case Title) TITLE filggk_1: (Control Data) TMIN, TMAX, NF, NCAL, MAXITR, IPNITR filggk_2: (Regularization Constants) OMEGAO, OMEGAl filggk_1: (Piecewise-Polynomial Specification) NSUB IDEG(1), NND(1) NDSUB(1,1), NDSUB(1,2), ..., NDSUB(1,NND(1)) NDBETA(1,1), NDBETA(1,2), ..., NDBETA(1,IDEG(1)+1) IDEG(NSUB), NND(NSUB) NDSUB(NSUB,1), NDSUB(NSUB,2), ..., NDSUB(NSUB,NND(NSUB)) NDBETA(NSUB,1), NDBETA(NSUB,2), ..., NDBETA(NSUB,IDEG(NSUB)+1) filggk_fl: (Filtering Information) IPRE, PREl, PRE2, PRE3 IPOST, POSTl, POST2, POST3 fllggk_§: (Sensor Information) NS XS(1), YS(l), W(1) XS(NS), YS(NS), W(NS) filggk_§: (Measurements) TIME(1), Y(l,l), ...., Y(NS,1) TIME(NT), Y(1,NT), ...., Y(NS,NT) 67 h.5 DISCUSSIONS ON SOME INPUT VARIABLES 0P INIHCP.DAI The purpose of this section is to provide background and insight into the assignment of certain input variables. Before discussing some of the input variables of INIHCP.DAT, three quantities, bias, variance, and mean squared error are discussed (pp. 153, Beck, Blackwell, and St. Clair, 1985). The bias is defined as the difference between the expected value of the estimated result and the true value of the parameter. It can be expressed as 0h, - '51.... - wk“)! (4.23) where ka is the bias caused by a biased estimator for the parameter. In the same equation, 5km is the true value, and fikm is the estimated A value of the parameter. The notation E(flkm) stands for the expected A value of the estimated result, fikm’ which also is the estimated value using exact temperature data. The variance of estimated value is defined as mm) - a {hekm - ammnz} (4.24) which compares the estimated value to the expected value of the estimated value. Had exact data been used, the estimated value is 68 equal to the expected value, and, as a result, the variance is zero. So, the variance is caused by the errors in the temperature measurements . The last quantity is the mean squared error which is defined as :2 - EH3 - a )2) ‘ (a 25) km km km ' It measures the difference between estimated and true values, and represents the combined effect of both bias and variance. The relationship between these quantities is {in ' Van.) + ”120,. (4.26) In general, the desired objective is to have a minimal value of the mean squared error. (in. 4.5.1 .HINIHAL.TIHE AND MAXIHAL.TIHE In some cases, only part of the measurements are used for the estimation of unknown heat flux. Program IHCP2D can skip the first part of measurement data, store some desired data, and ignore the rest of the measurement data. This option is controlled by the input variables TMIN, the minimum time, and TMAX, the maximum time. The variable TMIN specifies the starting measurement time that the temperature data will be stored in the program, and the variable TMAX 69 indicates when to stop the reading. The possible range for these two variables is 0 S TMIN < TMAX 5 cf (4.27) where tf is the maximum time for measurement data. 4.5.2 NUMBER.0P FUTURE TIME STEPS The symbol NF is the FORTRAN integer name for the number of future time steps, Nf, described in Section 4.2. According to equation (4.7), the measured temperatures at Nf future time steps are included in the estimation of parameters at present time. In other words, measurements, in, l s j 5 NS, from tm to tm+Nf-l have been used to estimate the parameters fikm' l S k s Np’ at tm. Assuming the measurement time step sizes, Atm's, are the same, increasing Nf will increase the bias and decrease the variance. That is, the estimated parameter histories will be smoother but probably more biased compared to the true values. The bias comes from the assumption of temporarily constant parameters during these Nf time steps, as described in equation (4.8). If the true curves for the parameters involve drastic temporal variations, this assumption can result in significant bias. On the other hand, if the temporal variation is gradual, negligible bias is expected. This feature has a similar effect as the first order regularization in time (Section 4.2). 70 This technique of using future Nf measurements can stabilize the algorithm as regularization does. The possible range for this number of future steps for B being held constant is 1 S N S 10 (4.28) in program IHCP2D. Experience has demonstrated that at least two future time steps are needed. Frequently used values are 2, 3, and 4. More are needed if the measurement errors increase in magnitude and also as the sensors are embedded deeper in the test specimen. The recommended number of future time steps is a function of the measurement errors and the sensor depth as discussed in Section 5.6 of Beck, Blackwell, and St. Clair (1985). 4.5.3 RATIO OF TIME STEP SIZES The ratio of time step sizes, NCAL, is defined as m: . NCAL - J z 1 (4.29) At c where Atc is the calculation time step size used by the direct problem solver to calculate temperatures. A proper value of Atc should be found to insure accuracy of the calculated temperatures. This can be done by trying several integer values (such as NCAL - l, 2, and 4) for a typical problem. A common value of NCAL is 2. 71 4.5.4 MAXIHUH.NUHBER.OF ITERATIONS The symbol MAXITR is the FORTRAN integer name for the maximum number of iterations allowed in finding the convergence of the parameter values. If this number is equal to zero, no iteration is allowed which corresponds to only one calculation. For linear problems (constant thermal properties and no radiation boundary condition), zero should be assigned to MAXITR to conserve computational effort. If the problem is highly nonlinear, l or 2 iterations may be needed to achieve convergence. In most nonlinear cases, MAXITR can be set to zero. If. there is uncertainty, the same problem can be solved using different values, such as l and 2, for MAXITR. 4.5.5 REGULARIZATION CONSTANTS The input variables OMEGAO and OMEGAl are the FORTRAN variable names for the zeroth and the first order regularization constants, mo and wl, discussed in Section 4.2. The zeroth order regularization has an effect of reducing the magnitude of parameters and damping the oscillatory behavior in the estimated values. This aspect can be seen from equation (4.16) if “1 - O and mo is assigned a large number compared to the entries in matrix A and vector C . Equation (4.16) then becomes I 50“”) - i7 0 (4.30a) 0 72 which indicates that 58*“ s o (4.301)) if ”0 approaches infinity. The first order regularization, on the other hand, tends to reduce the difference between subsequent parameters. This can be seen from equation (4.16) again, if mo - 0 and a large value has been assigned to wl. Equation (4.16) becomes “EH1 ”(M1) ' i: c (4.3la) 0r - ($119225... - 55.21%; .. which indicates that the differences in subsequent parameters approaches zero as ml approaches infinity. That is, it forces the parameters to be the same values. If these subsequent parameters are also spatially adjacent, this order of regularization causes the spatial variation of heat flux to be uniform. Selection of the values for mo and wl depends on the squared sum of residuals, which is defined as 73 N3 Nf 2 Squared Sum of Residuals - E) E (Yj,m+i-l - Tj,m+i-1) (4.32) j-l i-l Beck, Blackwell, and St. Clair (pp. 140, 1985) suggested that the regularization constant is chosen so that the squared sum of residuals is approximately equal to Nstaz, where a is the standard deviation of the measurements. In other words, the regularization constant is selected such that the root-mean-squares of the temperature residuals in the output file is approximately equal to the standard deviation of the measurements. Another formula for selecting this value is (pp. 211, Beck, Blackwell, and St. Clair, 1985) no or o1 a ("3“ wJ)(a/6,6)2 - (4.33) where (mix wj) is the maximum weighting constant of the sensors, 0 is the standard deviation of the measurement errors, and 65 is a measure of the change in the parameter from one time to another. However, a more direct approach is to try different values to study their effects on the estimated results. 74 Table 4.2 Values of Nf, u , and 01 for Stable Algorithm Nf “’o “’1 z1 20 20 >1 :0 20 21 #0 #0 75 In the case that the number of sensors, N8, is smaller than that of parameters, Np’ at least one of the regularizations should be activated to insure stability of the algorithm. Table 4.2 summarizes the values of Nf, mo, and ml for stable algorithm. 4.5.6 LOCATIONS OF SENSORS The FORTRAN variable names for the coordinates of the jth sensor, (x ,yJ), are XS(J) and YS(J), respectively. .The sensors should be J placed at locations which are most sensitive to the parameters. For the case that Ns - Np, each sensor should be placed at a location as close as possible to the parameter node. If the surface of the test specimen is in a harsh environment, such as high temperature or high flow speed, the sensors should be embedded slightly under the surface. For the case that Ns < Np' the sensors should be moved further into the solid so that the sensors can gather information from several parameters simultaneously (Hensel and Hills, 1987). A procedure for finding the most sensitive location involves calculation of the sensitivity coefficients. See Chapter 8 of Beck and Arnold (1977) for discussion of the optimal locations which can be found using these sensitivity coefficients. 4.5.7 UEIGHTING CONSTANTS OF SENSORS The corresponding FORTRAN variable name for the weighting constant of the jth sensor, w , is W(J). According to equation (4.7), the J 76 weighting constant represents the preference of the measurements taken at different sensors. For example, if all w - 0 but w1 fl 0, equation J (4.7) becomes Nf Np S - w (Y - T )2 + w 52 1 1,m+i-l 1,m+i-l 0 km 1-1 k-l N -1 p 2 + “1 E (pk+l,m ' 5km) (“'3“) k-l which indicates that only the measurements from the first sensor have been used for the estimation. A candidate for the value of the weighting constant is the reciprocal of the variance of measurement errors; that is wj - 2 (4.35) where aj is the standard deviation of the measurement errors in the jth sensor. It can be approximated by using the root-mean-squares of the temperature residuals. If the standard deviations of measurement errors are unknown, assign 1. to all weighting constants to bypass this feature. 77 4.6 EXAMPLE PROBLEM: TWO-DIMENSIONAL LEADING EDGE GEOHETRY In this section, program IHCPZD is employed to solve a two- dimensional problem. Pseudo-experimental temperatures are generated using a direct problem solver. These data then are entered to program IHCPZD to estimate the surface heat flux distribution. Three factors which affect the estimated results are studied. These factors are the number of future time steps, the number of iterations allowed, and a single error in the measurement. The geometry and nodal definitions of this example are shown in Figure 4.2. The spatial domain is denoted by 0. Part of its boundary is insulated and denoted by F1. The initial condition is the uniform temperature of To. Starting at a certain time, an unknown heat flux is applied to the boundary denoted by F2, and the temperatures of the solid begins to change. Radiative heat transfer on P2 also is considered because in practice the surface temperatures can be as high as 1500 K. The objective of this problem is to estimate the unknown heat flux applied to F2. Three sensors have been used to record the transient temperature measurements. If the origin is at node 2, the coordinates of these sensors are (0.0075 m,0.0035 m), (0.0015 m,0 m), and.(0.0075 m,-0.003S m), respectively. The sensor locations are shown and denoted by 31, 82, and 33 in Figure 4.2. The mathematical model of this problem is 511 1- (k.g§) + a— (k Q1) - pc 3: (x,y) in n, o s t s c 6x ay ay (4.36) f 78 Jenn -o (x,y) on r1, 05 cs: (4.37) f T(x,y,0) - T (x,y) in 0 (4.38) The material has constant density, and the thermal conductivity and the specific heat are temperature dependent. The material properties are given in Table 4.3. Property values of k1 and c1 between tabulated temperatures are calculated by linear interpolation. The transient temperature measurements are Y ji - T(x j'yj't1)+‘_11 15st,1515Nt (4.39) The measurement time step size is 0.1 s, and the number of measurement time steps is 30. The unknown surface heat flux q(x,y,t) on P4 is defined by _§_I -11 _ .4,“ k ax nx k 6y ny q(x,y,t) + as (T T”) (x,y) on F2, 0 s t S tf (4.40) where a is Stefan-Boltzmann constant, e' is the emissivity of the surface, and Tan is the ambient temperature. Notice that this is different from the usual inverse heat conduction problems because the not heat flux is not being found; rather the input heat flux is being found, assuming that the heat losses are only resulting from radiation. Also observe that the use of equation (4.40) makes the problem nonlinear. 79 1a 23 28 33 38 _ r” /7‘ P2 13 / X 8/ 51 / / F1 12 17 22 27 0 32 37 § E 6 11 / .4 2 —""' °. 32 16 21 26 31 36 / c 5 1o / 1 / \ 15 20 25 3o 35 / 4\ 3“ / 3 9 / \ JL 14 19 24 29 34 L *1 r. 0.01 m 4 _1 Figure 4.2 Geometry and.Nodal Definitions of Example Problem 80 Table 4.3 Material Properties and Some Constants Used in Example Problem T (K) 277.8 555.6 833.3 1111.1 1388.9 1666.7 1944.4 2222.2 c (J/kg-K) 397. 694. 865. 982. 1070. 1140. 1200. 1250. k (W/m-K) 185.2 148.8 122.0 103.8 90.5 83.4 79.6 77.9 p: 1862 kg/m3 6': 0.9 To: 294 x T : 294 x 0 R1 - 0.005 m R2 - 0.0157 m R3 - 0.005 m a1 - -o.25x107 w/n2 a2 - -1.x107 W/m2 a3 - -o.5x1o7 w/n2 81 The surface heat flux used to generate the simulated data is described below. First, the boundary P2 is divided into three subdomains, P 1 (from node 38 to 23), F2 2 (from node 23 to 19), and 2 P2 3 (from node 19 to 34). The heat flux on P2 1 is uniform. A second degree polynomial is assigned to the heat flux on F2 2. Another uniform heat flux is assigned to P2 3. The mathematical expression of this test case is q(r1,t) - 81 0 5 r1 5 R1, 0 s t s tf (4.41a) r2 r - - .2 - - .2 Q(r2at) 2(51 2fl2 + 33) R2 (331 “£2 + fi3) R2 + fil 2 0 5 r2 5 R2, 0 s t 5 cf (4.41b) q(r3,t) - 53 0 5 r3 5 R3, 0 s t s tf (4.41c) where 82 - q(‘2,t). The values of R's are listed in Table 4.3. According to Figure 4.2, these parameters, 81, 62, and 83, corresponding to heat flux components at nodes 23, 2, and 19, respectively. These parameters are time dependent, and their temporal variation are pk(t) - O 0 s s t s 0 5 s (4.428) fik(t) - ak (t - 0.5) 0.5 s < t S 1.5 s (4.42b) flk(t) - ak - ak (t - 1.5) 1.5 s < t S 2.5 s (4.42c) fik(t) - 0 2.5 s < t S 3 s (4.42d) 82 for k - l, 2, and 3. Equation (4.42) stands for a triangular test case. The constants a1, a2, and a3 associated with different fik's are also listed in Table 4.3. These constants are nggggixg so that the heat flux is flowing into the system. The temperatures of the solid rise from 294 K to 1400 R in 3 s, and the temperature difference within the solid is approximately 200 K during the peak period of the heat flux input. These calculated temperatures are used in program IHCP2D to estimate heat flux. First to be considered are the effects of the number of future time steps, Nf. According to equation (4.8), in the sequential method the parameters are assumed to be temporarily constant within these Nf time steps. However, the measurements used to estimate these constant parameters are generated from a test case with linear variation, equation (4.42). As a result, the estimated values can be interpreted as the weighted averages of the true values within these time steps. This difference between the estimated and the true values due to the assumption of temporary invariance is called the bias, which increases as the number of future time steps increases. The following presents an analysis of the bias due to different values of future time steps, which are Nf - l, 2, and 3. Exact temperature data are used in this study. Since the number of parameters is the same as sensors, regularization is not needed, and the results are obtained using the same piecewise-polynomial function as the test case. Figures 4.3, 4.4, and 4.5 show the estimated results of heat fluxes at node 23 (81), node 2 (£2), and node 19 (83), respectively. 83 In the case of Nf - 1, the results have very good agreements with the true curve because neither the measurement error or the bias is involved. This indicates that the inverse heat conduction algorithm is capable of producing accurate results. The curves of Nf - 2 and 3 also match the true values closely. However, in Figure 4.4 the curve of NE - 3 has approximately 7% bias (0.7x106 W/mz) in the region of linear variation. This is because the true values of three consecutive components vary 20% (0.2x107 W/mz) in this region. Using smaller time step sizes or the heat flux variation being smoother will ease this problem. In general, the results confirm that an increase in the number of future time steps increases the bias. The sign of the bias depends on the variation of the true curve. For example, if the heat flux increases in time, the estimated results will be greater than the true values and vice versa. Next topic of interest is the effects of the number of iterations allowed in the algorithm. Usually, several iterations are required to calculate the parameters for nonlinear problems (Figure 4.1). In this analysis, three cases are considered, in which the number of iterations are set to 0, 1, and no limit. (Although no limit is used in the third case, the parameters usually converge within 3 or 4 iterations.) Again, exact data are used in the estimation, and the number of future time steps is 2. The estimated heat flux at node 23 (pl), node 2 (52), and node 19 (83), are presented in Figures 4.6, 4.7, and 4.8, respectively. The results of these cases are indistinguishable, which implies that 84 iteration is not necessary for this problem in which the temperatures rise 1100 K in 3 s and a radiation boundary condition is involved. This aspect can be explained better by examining equation (4.15). In pm) T(n) , J m+i-1's’ and X(n)'s are this equation, the quantities, jik , updated during each iteration. If the results from 0 iteration are very close to those from n iterations, the values, 6(0), T§0;+1-1'S, and Xéga's, must be close to ”(n)’ T§?;+1-1's, and X§?;'s. In other words, the initial guess, 6(0), is near the estimated value. Based on this observation, if the change in the estimgggd_gglgg§ (not the true values) from time step to time step is small the number of iterations will not have a noticeable effect on the results. Last to be considered is the effects of measurement errors on the estimated results. This aspect is studied by introducing a single error of 20 K into measurement at 1 s in different sensors. Figures 4.9, 4.10, and 4.11 present the estimated results due to a single error in sensor 1, sensor 2, and sensor 3, respectively. In these cases, the number of future time steps is 2, which stabilizes the algorithm and is less sensitive to measurement errors. Figure 4.9 shows that the impulsive error in sensor 1 has most effect on the nearest parameter, #1, and is dispersed into one previous and two future time steps. However, its effect on two other parameters has been attenuated. This aspect can be explained using equation (4.15). In the right hand side of this equation, the residuals are weighted not only by the weighting constant wj but also by the sensitivity coefficient X112. Since parameter 81 is most close to sensor 1, the sensitivity coefficient 85 X111 is greater than the others. Hence, the estimation of fil weights heavily on the measurements from sensor 1. This behavior is also observed in Figures 4.10 and 4.11. Summarily, program IHCP2D has been used to solve a two-dimensional inverse heat conduction problem and demonstrated the ability of producing accurate results. It is found that (1) an increase in the number of future time steps increases the bias, (2) the number of iterations does not have a noticeable effect on the results, and (3) a single error in one of the sensor has the most effect on the nearest parameter. The CPU times for each of these cases are approximately 20 minutes on a VAX 750 computer. 86 MN ovoz um NSHH mac: woumluuam so macaw oluh oususm mo heels: mo nucouwm n.c ouzuuh any 08; 0..» new ohm mm.— CL _ 0.6 0.0 I oak-”LIMI- n Pb h n b n b n p p p u p n r t BO+UNo—.l r pnez o .. w .. I. Nnuuz D r: . ..l 3 .. mnez < .. No.50 _ Ma r r w r T m. m. Woo+modl W .. r H I r «u H. woo+moon .o. I - U an H .n. +8.1... x T r no N r! T / a - w n. loo+wo.~l two I r .l "a n T 87 N 3oz on 3E one: coon-Dam :o macaw mama. chasm mo Hon—la: .uc newsman so Gun—Mum 3 so: on mm ohm mm. om. mud od .1 . . Jaw Hr. - 1- . B+m~ T. H ......z o m 1 «1.2 a £38.... a. 1 ”"52 d o o r n...“ r . a T w .r 2 .. v . m. r ._ < . .... r mo+mo ml W I ... . r H I w ”- T a r . z 18651 m. r ... ._ u 1: .. . ... T m. .. ,, . T $3... x r u 1 / r ,. .. w .. m 18+mo.~1 b. a ., . a n n .. m. .... . _ p b . P _ . _ . . . L b . . m P: {POAV 88 2 one: on .33 23: noon-Sou :o macaw clan. chasm me Hi we nuoomum no cue-mum 3 2:: 9n mN 0N m; o.— ad od - h b r P P F L n — b n n p — n n p P h n n b h b n n b b ho+mmo F I... 1 03b. ....l T I Fnhz o r T Nnbz D r o ... rho+mo T. .. nuez o + . W ...- Woo+modl . T r T I Wmo+modl m r I Woo+moél m n r Tmo+uo.m1 T 1 1 ... .. r (55% (ZN/M) xnu 19914 99191141183 89 2 33: on .8: you: touwluumm :c accuuwuouu Ho koala! Ho wuoouum o.¢ Gunman 9n 9N bbrn—PhF-h-phPhbanP-P-P—D- A3 65.1.. ON m. — o. — m6 h ed on: onto: .15.. 835...: IIIIIIIIIIIIIIIIIIITII 16 C310 56.7w” P I ho+mo. Pl [Ulriri l00+modl Woo+uodl r r foo+mo.vl r Wmo+woNl (gm/M) xms 19914 9910141193 90 .4. so: on 89:. o8: comm-Jana :o accuuwuouu Ho hon—I52 Ho nucoumm h..~ 0H3»: A3 25» IL . p .035”. Llamas L . T NO+MN Fl .... OHCO: O W n .182. 8 (8+3... m .. nogEzc: 4 . T .... T . r w r . r O .... . ...-oo+mo.m1 w... T . . T H I. r. a T. . 18+moo1 m. T u H. T . T n T T +mo+ x T. r. I \I T . . T we M 1 .. r / 1 . 181mg: (a. T. . T. T _. .. T . .. s s e . his . . T. . . .1.Li_1t s . .glOAu 91 2 so: on as... one. tmumamuwm :o macaumuoun Ho HQAISZ Ho mucouum m.e 0H:&um 3 2:: o.n mmN OWN mm— ow— who od 1 p a, nmzflfiw pllHT. 4 L_ p p p p p b ~nh¥+mm~n F11 H ones. o M Tu. Fnohufi D IINTQTTMDoFII “.— .. 836:5 o T 4. . ... m m. fancied... n. F r H 1 1 .- TT ”T093061 w m T m m. .. + .T x - T 8 mo M w - m 1 ”Boreas... (w. e 92 m .— 1.. slug. mm a How—now a“ woman oumcum cu one ““051; was: woman—am cc wuoouum To 9.5»?— 3 2:: on 3 as m; 3 nd ed In PBML.~°L null-l.- — P P. -L n m n n — n n n b P n n n b Ih°+UNoFI n «12 .... H 1 no .82 o 18.3.5.7 a. 1 N ovoz D r ".1." w o. oooz o w m .... woo+mod1 w. r 1 H 1 1 G 1 18.58.? w H n u T T m T. r10 + . I1 \I . T o moo M T T. / T T w ..1 ”.8331 b. n h - lod 93 9n n a u oluh um N Howcom cu woman mamswm cu out moxsdm you: toad-Huwm no wuoommm On.¢ shaman gm 3 08:. ed n; o; no PinhLF-thbnbb—Pb-n—h-bnb- o.o TTTIUII'IIU'WIIITUIIIII tote oz «1.2 mm oooz a .82 o. 332 41:10]] samba. «.1 I to 18.58. T Woo+mod1 r- F Wco+ucdl r. Woo+uc.¢l Wmo+wodl (gm/M) x1113 10°11 99191111193 94 m A 1 club as n Manama cm Houum onwaam cu out moxzuh use: vouuluuwu cc wuooumu -.¢ Gunman 3 2:: on on o.~ n; o._ no PbLb—bLbrL—L-nnhbhbfibnbbb—h- 0.0 .95 oz «1.2 nu oooz a .82 o. 232 ‘°°ll lfrl—TFIITIIIFIIIIIIIIIU . T T BENT. r “.2137. r- .1oo+mo.m1 r- “T001000! moo+mo.¢1 1830a... (zW/M) an-1 1991-1 99191111163 - 1.0.0 95 4.7 OONCEUSIONS A computer program, IHCP2D, for the solution of general two- dimensional, nonlinear inverse heat conduction problems has been implemented and verified. It employs the combined method of function specification and regularization. In this particular method, a piecewise-polynomial function, which can be the piecewise combination of different degree polynomials, is used to approximate the spatial distribution of the unknown heat flux, and the function value at subdomain junction can be either discontinuous or continuous. The use of piecewise-polynomial function has greatly improved the flexibility of approximating the unknown heat flux. This combined method together with the piecewise-polynomial specification has never been reported before. An existing, powerful direct heat conduction problem code, TOPAZ, is modified into a subroutine and incorporated in the inverse heat conduction algorithm. The resulting computer program IHCP2D is capable of solving very general two-dimensional, nonlinear inverse heat conduction problems. This procedure of using the existing code TOPAZ gives greater power, generality, and reliability with less programming effort than if a two-dimensional direct problem solver were written specifically for inverse heat conduction problems. Although Hills (1987) has employed this concept to develop a parameter estimation code ESTIM, it is used for the first time to solve inverse heat conduction problems. At the present time, program IHCP2D is capable of solving two- dimensional inverse heat conduction problems with temperature dependent 96 properties, orthotropic materials, temperature boundary conditions, flux boundary conditions, convective boundary conditions, and radiation boundary conditions. More tests are required for problems with contact conductance, radiation in enclosures, and chemical reactions. Program IHCP2D has been used to solve a two-dimensional problem and demonstrated the ability of producing accurate results. The effects on the estimated results due to (l) the number of future time steps, (2) the number of iterations, and (3) a single error in the sensors are also studied. CHAPTER 5 SOLUTION OF.A ONE-DIMENSIONAL INVERSE MELTING PROBLEM 5.1 INTRODUCTION This chapter presents the solution method for the one-dimensional inverse melting problem. In this problem, both the temperature and the interface location measurements are used to estimate the unknown surface heat flux. A sequential function specification method is derived to utilize both types of information. The interface tracing method (Katz and Rubinsky, 1984, and Goodrich, 1978) is employed to solve the one-dimensional direct melting problem. A detailed procedure of this technique is described in Appendix B, and a direct problem code, MELT, has been written based on this method. Program IHCP2D has been modified by replacing SUBTPZ with MELT to solve the one- dimensional inverse melting problem. The following outlines each section in this chapter. The problem is described in Section 5.2. The sequential function specification 97 98 method, which is derived especially for this problem of different measurements, is presented in Section 5.3. The direct phase-change problem code, MELT, is introduced in Section 5.4. In Section 5.5, the sensitivity coefficients of the problem are investigated. The estimated results for the surface heat flux are compared to true values in Section 5.6. The effects on the estimated results due to (a) the number of future measurements used, (b) measurement errors, (c) pre- filtering and post-filtering, (d) experimental time step sizes, and (e) types of measurements are also explored. The last section concludes this chapter. 5.2 PROBLEM DESCRIPTION In this section, a mathematical description of the one-dimensional inverse melting problem is given. The geometry of the problem is shown in Figure 5.1. In the initial condition, the material is in a solid phase with a uniform melting temperature. An unknown heat flux is then applied to the left boundary of this one-dimensional finite slab. The right boundary is insulated, and the material begins to melt. The temperature measurements are taken at x - x8, and the propagation of the interface is assumed known. The objective of this problem is to estimate the unknown surface heat flux applied to the left boundary. The mathematical formulation of this problem is og—g-gi 0=H Hmcouucolu010cc Ho huuoloou ~.n unawuh \\\\ \\\\\\\ oootoE. 020m 233.. k. 102 The inclusion of the interface location history in the estimation procedure is significant. This means that different types of measurements quantities are used together to estimate the unknown function. Using this idea, other types of measurement can also contribute to the estimation procedure. One example might be the measurements of an interior heat flux. Another important point is that the inclusion of the interface location in the estimation procedure introduces the possibility of an inverse design problem in addition to an estimation problem. In some heat treatment problems, it is desired to specify the movement of a certain temperature front which is accomplished by prescribing an appropriate surface heat flux history. The analysis includes this design possibility. 5.3 SEQUENTIAL FUNCTION SPECIFICATION METHOD The sequential function specification method is employed to solve this one-dimensional inverse melting problem. The procedure derived below is similar to that in Section 4.2 except for the objective function. The unknown function q(t) is first discretized into Nt unknown parameters; that is q(t) - 8m tm-l < c s cm, 1 s 1 s N (5.11) I: These parameters then are estimated sequentially. The objective function for this problem is 103 N f 2 2 s ' E [ wTatum-1 ' Ts,Tn.+1-1) + "u(Yu,m+1-1 ' ”n+1-1) ] (5'12) i-l where Ts,m+i-l and "m+i-1 are functions of 5m' The set of weighting constants is defined with "T being the reciprocal of the variance of e and wu being the reciprocal of the variance of 6 Ti vi' To estimate fim’ the Gauss method of minimization is employed to minimize the objective function. Following the same procedure described in Section 4.2, the resulting equation is fi§n+1) _ pén) + (n) (n) (n) (n) VT} xT,M “l( T, m+i- 1 Is, m+i- l) + “9} xv,m+i-1(Yv,m+i- lum+i-l) (0) )2+ n) W} (x1, m+i 1) R} U4 m+i 1)2 (5.13) where the superscript 'n' is the iteration number and the summation is (n) ”(n) (n) for i from 1 to Nf. In this equation, T8 1 1, "m+i-1’ xT,m+i-1’ and X(n) all are evaluated at fi(n) and the sensitivit ffi i I1-1 m , y coe c ents, x'r,n+1-1 a“d xv,m+i-l’ are 81 __.§...Ei1.'_l x1 m+i 1 afim 1 s 1 5 NE (5.14) By - + - xv,m+i-l 532—1—1 1 S i S Nf (5.15) 104 Comparing equations (5.13) and (4.15), it is found that the inverse heat conduction algorithm for the solution of two-dimensional IHCPs can also be used to solve this problem if both regularization constants are set equal to zero and the weighting constants have been properly selected. The standard deviation of heat flux, oq, can be approximated by 2 2 2 2 2 2 2 wTOT E xT,m+i-l + wuav E xu,m+i-l a g (5.16) q (w 2 + w X2 )2 T xT,m+i-1 u,m+i-1 (pp. 248, Beck and Arnold, 1977) where the summation is for i - 1 to Nf. The standard deviation of the temperature measurement errors is denoted by 0T, and the standard deviation of the interface measurement errors is denoted by 0". Equation (5.16) states that the variance in the estimated results is a weighted combination of the variances in the measurements . 5.4 IMPLEMENTATION OF MEET TO PROGRAM IHCPZD A computer program, MELT, has been developed to solve phase change problems. The interface tracing method (Katz and Rubinsky, 1984, and Goodrich, 1978) is employed in this code so that the interface location can be calculated explicitly. A detailed derivation of this method is presented in Appendix B. 105 Program IHCP2D is modified to solve this inverse melting problem. The major modification is to replace SUBTPZ, the direct problem solver for two-dimensional heat conduction problems, with MELT. The resulting computer program then is capable of solving the inverse problem described in Section 5.2. 5.5 SENSITIVITY COEFFICIENTS The sensitivity coefficients for this problem have been defined in equations (5.14) and (5.15). From the definition, a sensitivity coefficient is proportional to the change in T (or ”m+i-l) due s,m+i-l to a change in the value of fim' Since fin is assigned to q(t) over the range of tm-l < t s t this definition can be interpreted as the m+Nf-1’ change in T8 m+1_1 (or "n+1-1) due to a 9282.2hanss in q(t) fr°m tm-l to tm+N _1. This aspect of step change can be seen in the following f equation TE m!i-](fim+5fi) - '1'E m!I-](pm) x'r,n+1-1 “ 58 ' TR m'i-l(q19~9ql_lnqm+6fivfi|1+6:é-o) ‘ TW1-1(911-Tqm19m+11u) (517) Hence, this quantity actually is a step change sensitivity coefficient. A symbol tsc’ which will be used later, denotes the starting time when the step change occurs in q. In equation (5.17), the step change 106 starts from qm, so tsc is equal to tm-l because qm is assigned to heat flux during tm-l < t s tm. If the sensitivity coefficients are zero, T and u are s,m+i-1 m+i-l independent of fin, and have no value for estimating 5m. Thus, the study of the sensitivity coefficients is important to the design of the experiment. It suggests which, where, and when the measurements should be taken in order to effectively estimate the parameters. In order to illustrate the concept of the sensitivity coefficient, a particular example is considered; the heat flux function in this example is q(t) - l O s S t S 0.15 s (5.18a) q(t) - 1 + 135—2-43- o.15 s s c s 0.75 s (5.18b) q(t) - 2 0.75 s s t (5.18c) The dimension of q is W/mz. Figures 5.2, and 5.3 show the sensitivity coefficient histories of v and T m+i-l with respect to flm for s,m+i-1 different starting times tsc - 0.12 s, 0.42 s, and 0.72 s. For these cases, the sensitivity coefficients during tsc 5 t S tsc+0.3 s are evaluated. Figure 5.2 shows the sensitivity coefficient histories of u with respect to 5m; the curve of tsc - 0.72 s is relatively smaller than the other two curves. This is because at later times the interface moves 107 farther away from the surface x - 0. As a result, the measurements of interface locations become insensitive to the surface heat flux. In Figure 5.3 which shows the sensitivity coefficient histories of T(xs,t) with respect to fim’ it is observed that the curve of tsc - 0.12 s does not respond until t - 0.21 s. This is because that the interface did not reach x8 (- 0.2 m) until this time. Note that in this problem the initial condition is a uniform distribution at the melting temperature. As a result, a sensor located at xs - 0.2 m is not able to detect any activity to its left before the interface passes xs. If the sensitivity coefficient is a function of the parameter, the sensitivity coefficient history will depend on tsc (pp. 22-30, Beck et a1., 1985). Figures 5.2 and 5.3 reveal that both sensitivity coefficients are functions of the parameter pm. This type of estimation problem is considered nonlinear with respect to 3m“ 108 a mo monsoon; 0583030 #53383 «a 3:»: Amy Etc. «4 0.0 0.0 I . . _ . .TT....T..T!..T 1 f. 1 T T a I I 1 ‘L I 1 D l. 1 D II E I VL 1 a 1 a ..----- .0.- n. 0.0 ......e I C I I 2 I ‘\ I ‘5 I l a r. I C ‘5 1 F - T. «gnome To .. T H m meduoma mlm T . on m n... L 0.6 T T Tmmo_ .T _ T t. _ . . . _ . . . 00.0 (M/gw) 89/49 109 385... no 83883: asoaoaumooo 303382." an ohms.— 3 9:: N; no so no L b If _ h P I 1. ... I .. I 1: .1 1 .. .1. .l. .1. .1. . 1 .I. .1. 11111111 00.0 I- . = C T. . 2 c T. . r e a II . 2 e E O 1 o a 1 . a 1 a - o «honour To .. a $6qu 818 .. a «gnome I P p r F L . p _ b _ _ — . _ .l 00.0 19 (M/zw-x) 99/ 110 5.6 RESULTS AND DISCUSSIONS The surface heat flux described by equation (5.18) is used to test the inverse melting code. Pseudo-experimental data are generated using the direct problem code, MELT. These data then are entered in the inverse problem code to estimate the surface heat flux. Figure 5.4 shows the comparison of the estimated results (using errorless data) with true values. They display very close agreement for both cases of Nf - 2 and 4. There are two types of errors which occur in the estimated values. One is due to the bias caused by the use of future measurements. An average measure of the bias, D, is N-N+1 T. t f 2 2 ”-1 [fin ' 3(pm)] D - (5.19) Nt-Nf+1 A where E(fim) is the estimated value resulting from errorless data, and 5m is the true value. The other error is caused by the measurement errors. An estimate of the standard deviation in heat flux is N-N+1 T. cf "2 m-l Ifim ' E(fl&] a _ - (5.20) q Nt Nf + l ) N 111 A where fin denotes the estimated value resulting from data with random errors . The effects on the estimated values due to the following factors are of interest in this investigation. They are the number of future measurements Nf, experimental errors, pro-filtering (the measurements) and post-filtering (the estimated results), time step sizes, and the types Of measurements. First considered are the effects on the estimated results due to the number of future measurements Nf. In other words, it is the bias caused by Nf to be investigated. Four cases are considered, which are Nf - 1, 2, 3, and 4. From Table 5.1(a), it is seen that the bias increases rapidly as Nf increases. This suggests not using value of Nf greater than 4 or 5. Second considered are the effects on the estimated values due to the measurement errors. A random number generator has been used to super-impose additive errors on the calculated data to create pseudo- experimental data. The errors are normally distributed with constant standard deviation, 0", where n denotes either T or u. Three cases are considered, which are an - 0'01nmax’ 0.02nmax, and 0'03nmax' They are compared to the estimated results using errorless data. The symbol "max represents either I , the maximal temperature at x , or u , the s,max s max maximal interface location. Table 5.1(b) shows that the estimated 112 A standard deviations, aq's, corresponding to these three cases are linearly proportional to the standard deviation of these small errors. This can be explained by examining equation (5.16). It shows that when a and av change with the same proportion, aq also changes with the T same proportion. Third considered are the effects resulting from filtering processes. For pro-filtering process, the measurements are filtered using I "1 (5.21) - 0.25Y + 0.5Y + 0.25Y ’7 '7 '1 ,i-l i ,i+1 where n denotes either I or v. The post-filtering process filters the estimated parameter values using 3m - 0.258“.1 + 0.58m + 0.258n+1 (5.22) In equations (5.21) and (5.22), the symbol '~' denotes the smoothed quantities. Table 5.1(c) shows that introducing of either filtering A process reduce aq approximately by 40%. It seems as though the filtering process can effectively reduce the variance in the estimated data. In fact, this is because for the case considered the frequency of the measurement errors is the same as the sampling rate. In other words, there is no correlation in time in the the measurement errors. 113 Fourth considered are the effects of the measurement time step sizes, Atm. Four cases are considered, which are Atm - 0.02 s, 0.03 s, A 0.04 s, and 0.05 s. The estimated standard deviations, aq's, corresponding to these cases are shown in Table 5.1(d). The results do A not display a clear trend in the values of aq while increasing Atm. This can be explained as below. Figures 5.2 and 5.3 show that the step change sensitivity coefficients increase monotonically. This suggests that larger AtIn values result in larger sensitivity coefficient values. From equation (5.16), this increase in the sensitivity coefficient decreases the standard deviation aq. On the other hand for the same va larger Atm values also causes more bias due to the temporarily constant assumption of heat flux, equation (4.8). These two contrasting factors help to explain the results in Table 5.1(d). The last factor to be considered is the effect resulting from the type of measurements employed. The heat flux is estimated using the interface measurements and temperature measurements separately. According to equation (5.12), this can be accomplished by setting the weighting constants, "T and w”, to 0 and a;2 in one case and 0&2 and 0 in another case. The procedure fails to estimate the heat flux by using the temperature measurements alone, because the sensitivity coefficients are zero before interface passes through location x8. From Table 5.1(e), the interface measurements also result in large variance in the estimated values. If the variance in the interface measurements is decreased, it then becomes possible to estimate the heat flux using the interface data alone. 114 00.0 bblPLbLb—p—nPbb—nhprbbbP moa~w> mama saw: wuuawom vouwluunm Ho flomuumnlou ¢.n ousm~m 050 00.0 3 2:: $10 00.0 mp... 0 .00.0 I I III‘ITIIIIIIIIIITII on: I Nakz O 41.2 a ‘I C (r Iv Ir Ir II III Q I II‘ IO‘ pnb-_-—bP—I——F——L-b__PF-_Lpn— mud (zw/M) (1)b an 10:31.1 aooyns 115 TABLE 5.1 Summary of the Investigation.Resu1ts (a) W (not filtered, Atm - 0.03 s, both measurements used, exact data) Nf l 2 3 4 Bias (W/mz) 0.00002 0.00567 0.0126 0.0215 (b) Exnsrimental_firrsrs_Analxsis (Nf - 3, not filtered, Atm - 0.03 s, both measurements used) 0' -41 0.01 0.02 0.03 "m aq (W/mz) 0.0588 0.118 0.176 (c) Eilss:ins.£rgse§ses_Analxsis (Nf - 3, a" - 0'01nmax’ Atm - 0.03 3, both measurements used) Filtering No Pre-filter Post-filter aq (W/mz) 0.0588 0.0346 0.0384 116 TABLE 5 . 1 (Continued) (d) neaaurement_Iime_§ten_iizes_Analxsia (Nf - 3, an - 0.01nnax, not filtered, both measurements used) Atm (s) 0.02 0.03 0.04 0.05 aq (W/mz) 0.0925 0.0588 0.061 0.0307 (e) Measurement.1xnes.énal¥si§ (Nf - 3, a" - 0.01nmax, not filtered, AtIn - 0.03 3) Type Both Yui aq (W/mz) 0.0588 0.778 117 5.7 SUHHARY'AND CONCLUSIONS The sequential function specification method has been employed to solve the inverse melting problem. Although constant thermal properties were used, the procedure can handle temperature dependent properties as well. Both the temperature and the interface measurements are used to estimate the surface heat flux. This signifies that different types of measurements can be used together to improve the estimated results. Program IHCPZD has been modified to solve the inverse melting problem and demonstrated the ability of producing accurate results. The effects on the estimated results due to a) the number of future measurements used, b) experimental errors, c) pre-filtering and post-filtering processes, d) measurement time step sizes, and e) types of measurements are investigated. This study concludes that (1) increasing the number of future measurements improves stability, reduces the deleterious effects of measurement errors, but increases the bias; (2) the estimated standard deviation of heat flux are linearly proportional to the standard deviation of the measurement errors (for small errors); (3) both pre-filtering and post-filtering processes effectively reduce the estimated standard deviation of heat flux and have similar effect on the results; (4) increasing of the measurement time step size reduces the standard deviation of the estimated values, but increases the bias; and (S) for this particular problem, it is preferable that both measurements of temperature and interface location are used. CHAPTER 6 CONCLUSIONS AND FUTURE STUDIES This chapter concludes this dissertation, and some future studies are suggested. 6.1 CONCDUSIONS The function specification method for the solution of inverse heat conduction problems has been discussed in detail. It is found that numerous techniques can evolve from this method. Each is different from others in the approaches of approximating the unknown function, estimating the parameters, defining the objective function, minimizing the objective function, and calculation of the sensitivity coefficient (or gradient). These different combinations result in a variety of function specification methods. A detailed discussion of this method as presented in Chapter 2 has never been reported before. 118 119 It is also found that a main concept of the conventional function specification method concerns approximating the unknown heat flux, while the essence of the regularization method is to incorporate smoothing term(s) in the objective function. So, it is possible to combine these two widely used methods into one procedure. Three approaches are described for the calculation of the sensitivity coefficient, aT/afl, or gradient of the objective function, aS/afl. The first one is the direct method in which the sensitivity coefficient (or the gradient) is calculated by solving the sensitivity problem. The second approach is the adjoint method in which the sensitivity coefficient (or the gradient) is calculated indirectly through the adjoint variable. The third approach is the finite difference approximation in which the partial derivative is approximated by finite difference equation. For nonlinear inverse problem of estimating Np parameters using NS sensors, the direct method requires the solutions of Np linear sensitivity problems plus one nonlinear temperature problem. The adjoint method requires the solutions of N8 adjoint problems plus one nonlinear temperature problem. The finite difference method requires the solutions of Np+l nonlinear temperature problems. The choice of method relies on the relative magnitudes of Np and N3. If Np > N5, the adjoint method is preferable. If Np 5 N8, the direct method is preferable. However, extra effort in programming might be required because the sensitivity problem is different from the temperature 120 problem. For linear inverse problem with Np s Ns’ both the direct method and finite difference method can be used. Three different adjoint methods have been presented in this study. Two of them are for the calculation of the gradient and the other one is for the calculation of the sensitivity coefficient. Each has defined the adjoint problem differently to suit its own purpose. There has never been a publication that compares different adjoint methods as presented in Chapter 3. The first problem of interest in this dissertation is the solution of a general two-dimensional inverse heat conduction problem. A combined method of function specification and regularization has been employed for the solution of this problem. This method uses the piecewise-polynomial function, which can be any combination of different degree polynomials, to approximate the spatial distribution of the unknown heat flux. This piecewise-polynomial specification has greatly improved the flexibility of assuming a functional form for the unknown heat flux. Moreover, possible discontinuity in the unknown function can also be accounted for using this approach. This combined method together with the piecewise-polynomial specification has never been reported before. The direct problem solver of this program is adapted from an existing general two-dimensional heat conduction code, TOPAZ. As a result, the program is capable of solving two-dimensional nonlinear inverse heat conduction problems with great diversities. This gives greater power, generality, and reliability for the same effort than if 121 a two-dimensional direct code were written specifically for the inverse problem. This concept of using an existing direct problem code is for the first time employed to solve inverse heat conduction problems. A similar procedure can be applied to other types of engineering problems for function or parameter estimation using proper direct problem solvers. Program IMCPZD has been used successfully to solve a two- dimensional inverse problem. The effects on the estimated heat fluxes due to several factors are studied. It is found that (1) an increase in the number of future time steps increases the bias, (2) the number of iterations does not have noticeable effect on the results, (3) a single error in one of the sensor has most effect on the nearest parameter. The second problem of interest in this study is a one-dimensional inverse melting problem. A computer program, MELT, based on the interface tracing method has been written to solve the direct melting problem. By replacing TOPAZ with MELT, program IHCPZD is used to solve the inverse melting problem. Although constant thermal properties has been assumed in the problem description, the procedures can handle temperature dependent properties. Both the temperature and the interface measurements have been used to estimate the surface heat flux. This signifies that different types of measurements can be used together to improve the estimated results. This inclusion of different type of measurements is an unique characteristics of this study. This one-dimensional inverse melting code has demonstrated the ability to produce accurate results. The effects on the estimated results due to a) the number of future measurements used, b) experimental errors, c) pre-filtering and post-filtering processes, d) 122 measurement time step sizes, and e) types of measurements have also been investigated. 6.2 FUTURE STUDIES In this section, several possible studies in the area of inverse heat conduction problem are presented. One topic of interest is the development of a computer code for the solution of general three-dimensional inverse heat conduction problems. One major difference between a two-dimensional and a three- dimensional inverse problem codes is the direct problem solver. A three-dimensional direct problem solver should be used in the three- dimensional inverse problem code. Besides this, the functional form to approximate the unknown heat flux is also different. -In two- dimensional case, the unknown heat flux is a function of one surface coordinate, r. For three-dimensional problem, two spatial coordinates are required to accomplish this task, and the surface approximating functions, such as bicubic spline, should be used instead of piecewise- polynomial functions. Another possible project is to study the optimal locations of sensors, (xj,yj), the optimal number of sensors, Ns’ the optimal measurement time step sizes, Atm, and the optimal number of future measurements used, Nf. This problem can be posed as an optimal design problem. The objective function for this optimization problem is a norm calculated from the matrix of sensitivity coefficients, and the design variables include Ns’ x ( 1 s j S Ns)’ Atm, and N . 1' ’1' f 123 Numerical analysis can also be employed to study the errors, stability, convergence, and consistency of inverse problem algorithms. The results from this systematic study will be very helpful in selecting the inverse problem procedure. APPENDICES APPENDIX.A SUBROUTINE HIEEARCHY This appendix discusses the subroutine hierarchy (or program structure) of program IHCPZD. Figure A.l shows the subroutine hierarchy of program IHCPZD. The main program is IHCP2D. It first calls subroutine INPUT9 to read the input files. Given the nodal numbers of parameters for each subdomain, the subroutine INDEX then re-numbers the double subscript parameters, 511(t), to become single subscript parameter fik(t) and calculate the total number of the parameters, Np. Next, the main program calls COORD to construct a local coordinate system, rj, for each subdomain as indicated by equation (4.1). In code IHCP2D, a sensor location is not restricted to coincide with a nodal point. This feature allows mesh generation programs to be used for automatic mesh generation. The subroutine FINDND is used to search for the elements containing the sensors, then the sensor temperatures are interpolated from these element nodes. 124 IHCPZD OORD F FINDND F 4 LJ PREFIL ‘-~FILTER 'SUBTPZ ' NONLI [ OUTPU9‘ Fl [‘1 F! El I board-I FILTER —. J { ‘F""“ r-*PRHEAT 125 ‘ CALIEMT' SENSIT F0 GAUSS PRBETA CALTEM 'SUBTPZ UPDAT9 SUBTPZ TS 'DETERM3 _I nu DETERM3 DETERMA Figure A.1 Subroutine Hierarchy of Program IHCPZD 69:1 DETERM3 126 If pre-filtering is requested, the subroutine PREFIL is used to smooth the measurements. Next, the main program uses subroutine NONLI, the inverse heat conduction algorithm, to estimate the parameters simultaneously in space and sequentially in time. If post-filtering is requested, the subroutine POSTFL is used to smooth the estimated results. The last task is to print out the results which is accomplished by the subroutine OUTPU9. The subroutines of the second level are discussed in the following. The subroutine INPUT9 calls the subroutine SUBTPZ to read the input file INTPZ.DAT. The inverse heat conduction algorithm, NONLI, calls five other subroutines to estimate the parameters. The first one is CALTEM which calculates temperatures given guessed parameter values. The sensitivity coefficients are calculated-using SENSIT. Given these quantities, the system of equations is constructed using FORMAB. These equations are solved by subroutine GAUSS which utilizes Gauss-Jordan elimination method (pp. 272, Carnahan et a1., 1969). Subroutines PREFIL and POSTFL both call FILTER to perform the filtering process. Subroutine PRHEAT is used by OUTPU9 to print estimated nodal heat flux. The following describes the subroutines of the third level. To calculate the temperatures, subroutine CALTEM calls UPDAT9 to update the initial conditions and heat flux in the associated common blocks. Then, subroutine SUBTPZ is used to calculate the temperature distribution. With this information, the sensor temperatures are interpolated from the neighboring nodal temperatures by function TS. Function FCTQ is used by UPDAT9 to calculate nodal heat flux according to the Lagrange formula, (equations (4.4) and (4.5)). Functions 127 DETERM3 and DETERMa are used by TS to calculate the determinants of 3x3 and 4x4 matrices, respectively. APPENDIX.B INTEREACE TRACING METHOD FUR.TUE SOLUTION 0P ONE-DIMENSIONAL.PUASE CHANCE PROBLEM 8.1 INTRODUCTION Several methods have been reported on the solution of phase change problems including the analytical method (pp. 406, Ozisik, 1979), integral method (pp. 416, Ozisik, 1979), moving heat source method (pp. 423, Ozisik, 1979), enthalpy method (Voller and Cross. 1981), heat capacity method (Bonacina et a1., 1973), and interface tracing method (Katz and Rubinsky, 1984, and Goodrich, 1978). The interface tracing method will be used in this appendix to solve one-dimensional phase change problems. The following describes the outline of this appendix. The mathematical formulation of the problem is presented in Appendix B.2. The solution technique, interface tracing method, is described in Appendix 3.3. A detail derivation of the finite difference equation for the interface node is depicted in Appendix B.4. Appendix 8.5 presents the comparison of numerical and analytical results. 128 129 loanoum mcuudoz decedenolua1066 no hwuoloou ~.n ouswum \\\\\-\\\\\\\ a . _ _ fl _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ u _ n n _ _ _ _ _ _ I 3 . _ . _ _ _ _ _ _ oootofi: _ .1 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ u _ _ _ . . , _ _ . I _ . u _ . _ _ + z _ :6: _n_ _ _ m _ _ _ _ _ _ _ _ . _ _ _ n u _ _ _ _ 26m ” n _ _ 233.. u _ _ _ _ _ _ _ _ _ _ _ _ u _ _ _ . m _ . _ u n _ _ 3c 130 8.2 PROBLEM DESCRIPTION Figure 8.1 shows the geometry of a finite slab. The slab initially is at solid phase and the initial temperature distribution, f(x), is less than melting temperature which is assumed to be zero. Starting at time zero, a heat flux, q(t), is applied to the left surface, and the right surface is insulated. The purpose of this direct problem is to calculate the temperature distribution as well as the interface propagation. The mathematical formulation of the problem is 8T 3T L _ _& 3x (k1 8x ) pct at 0 < x < u(t), 0 s t 5 cf (8.1) a. 33: a}; ax (ks ax ) - pcs at u(t) < x < a, O s t 5 cf (B.2) Ts(x,0) - f(x) < O O s x s a (8.3) 3T - k1 a—xl (0,t) - q(t) 0 s c 5 cf (B.4) 8:. ax (a.t) - 0 O s t 5 cf (8.5) The subscripts 's' and '1' stand for solid phase and liquid phase, respectively. According to equation (B.4), the heat flux is positive if entering the slab. The interface conditions are T2(V,t) - Ts(v,t) - O O S t S t (8.6) f 1 ax 3 3x ' PL dt x - u(t), 0 s t s tf (8.7) 131 It is assumed that there exists a discrete line of interface and the melting also occurs at a distinct temperature rather than a region of temperature. The thermal conductivities, ks and k2, and the specific heats, c8 and c can be function of position. The symbol L stands for l! latent heat which is a constant. The densities of both phases are assumed to be the same to neglect the mass convective effect. 8.3 INTERFACE TRACING METHOD The interface tracing method used to solve the problem described by the equations (8.l)-(8.7) is described below. First, the domain is divided into N equally spaced grids with one node located on each grid line (Figure 8.1). For convenience, the temperature for node j at time tm is represented as T(x ,tm) - T (8.8) J 1!!! Next, a control volume is assigned to each node, and the control surface is defined half way between two adjacent nodes. One extra node is defined at the interface location and travels with it. This moving node is called the interface node and designated by 'f' in Figure 8.1. Its adjacent stationary nodes are designated by 'p' and 'p+l'. The control boundaries between node f and p to the left and p+1 to the right are also defined as half way between them. If the interface node is moving from left to right with the speed du/dt, these control boundaries are also moving but with half of the speed. As a result, the control volumes of nodes f, p, and p+l are deforming. In 132 other words, an interface node with finite but deforming control volume has been created to trace the interface location. The energy balance equation for a deforming control volume (pp. 121, Potter and Foss, 1975) is g?! pudV-I-I puv-ndA-i-I E-BdA-o (8.9) CV CS CS where the subscripts 'cv' and 'cs' stand for control volume and control surface, respectively. The density, internal energy, and the relative velocity of mass flux to the control surface, are represented by p, u, and 3, respectively. The volume and area are denoted by V and A, respectively. Heat flux applied to the control surface is represented by 3, and a is an unit vector whose direction is outwardly normal to the control surface. The first term on the left of the equation is the rate of change of internal energy within the control volume including the effect of moving boundary. The second term accounts for the internal energy carried out of the control surface by the mass flux. The third term accounts for the total heat flux out of the control volume through the control surface. This energy equation for deforming control volume is then applied to the control volumes of nodes p, f, and p+1. The energy equations for these three nodes can be solved for three unknowns, T , T pm p+l,m’ and um. The same energy equation can also be applied to the control 133 volumes of other stationary nodes. A detailed derivation of the finite difference equation for the interface node is presented in Appendix 8.4, and the equations for all the nodes are summarized in Table 8.1. With calculated rpm and T backward substitution of p+l,m’ elimination algorithm for tridiagonal system (Pp. 154, Conte and de Boor, 1980) is applied to the nodal equations from nodes 0 to p-l to The forward substitution solve for the temperatures T T p-l,m .... Om' is employed to the nodal equations from nodes p+2 to N to solve for the temperatures Tp+2,m .... TNm' 8.4 DERIVATION OF FINITE DIFFERENCE EQUATION FOR INTEREACE NODE This section derives the finite difference equation for the interface node. For convenience, backward difference approximation is employed to approximate temporal derivatives. The control volume of the interface node is shown in Figure 8.1. Applying equation (8.9) to this control volume, the first term on the right hand side becomes (u+x )/2 t dt cv (V+Xp)/2 where the volume integral has become a linear integral. For a melting problem, the internal energy of liquid phase also includes the latent heat, L. Since the interface is moving toward the right, this control volume is filling up with the latent heat just like a bucket. The internal energy is divided into two terms. One is latent heat, and the other is the specific heat as 134 (v+x )/2 LI pudV-A[LI 9+1 chdx+§E pde] (8.11) cv (v+xp)/2 (v+xp)/2 To perform time derivative on a deformable control volume, Leibnitz's rule is employed which is (t) g; F(X:t) dx - ah F(b; t) + J'g— F(x; t) dx - g5 F(a; t) (8.12) a(t) The equation (8.11) then becomes v+x v+X g; pu dV - A [ % fig ch(-§nil,t) - 2 1%? ch(-—n, t) cv (v+x )/2 ° + p+l- L ch dx + Q PL - A in (8.13) (v+x )/2 at dt 2 dt P The second term on the right hand side of the equation (8.9) is +x v+x *.* - 112 :__n - lgg + LN“ [... w] ...»cu—Ll and the third term is +3 u+x ...... - flr—n _ n n+1 Icsq n dA A [ k ax( 2 ,t) k ax( 2 ,t) ] (8.15) Substitute these equations into (8.9), it becomes (u+x )/2 u+x v+x P” pcfldx+1¥pL+kfl(—Jt)-kfl(—Eflt)-o (8.16) 9 Finite difference approximation is used to approximate the partial derivatives in this equation, and it becomes pcp(vn+1 - xp) (Tp,m+1 - Tpm) + pcp+1(xp+1 - Vm+1) (Tp+l,m - Tp+1.m) k 'r k '1' + + + 8pL(Vm+1 - um) - saga-1313;- + W _ y ) - 0 (8.17) m+1 p p+l m+1 where op represents the specific heat at the pth node, and the melting temperature has been assumed to be zero. This is the finite difference equation for the interface node, and it allows position-dependent properties. 136 TABLE 8.1 Finite Difference Equations for all the Nodes Nede_fi pcp(um+1 - Xp) (Tp,m+1 ' Tpm) + P¢p+1(xp+1 ' "m+l) (Tp+l,m - Tp+1.m) k T k T ' _ _ _n_nim:l_ .n:l_n:lim¢l _ + 8pL(um+1 um) 88c(y . x + x _ y ) 0 (8.18) m+1 p p+l m+1 8242.2: 3("m+1 - xp-l) (Tp,m+l ' Tpm) + (xp . xp-l) (Tp-l,m+l Tp-l,m) 39.8181 ‘ TD-lmfl MT + 8 Ata ( x _ x ‘ + y - x ) - 0 (8.19) p p p-l m+1 p If 2 - Q: T .2411 211:1 3ym+l(T0,m+l - TOm) + 8Atao( ”m+1 - k0 ) - 0 (8.20) uggg n+1: 3("p+2 ' ”m+1) (Tp+l,m+l ' Tp+1,m) + (xp+2 ' “p+1) (Tp+2,m+l ' TP+2,m) r r - r p+l,m11 p+2,m¢1 n+1,m¢1 _ + 8Atap+1(x _ u + x - x ) 0 (8.21) p+1 m+1 p+2 p+1 137 TABLE 8.1 (Continued) 8 ..Nimil__ “"19 ' ”m+1) (TN,m+l ' TNm) + ““19 xN - um“ " ° Ng§g_Q: (not nearby the interface) a At a At a At T0 m+1(1 + 2 2) ‘ 2T1 m+1 2 ' T0 m+1 + 2 2 qm+1 ’ Ax ' Ax ’ Ax I121§a1_flggg_1: (not nearby the interface) T1-1 m+1 a A: + T1 1181‘1 + 29 A2) ' T1+1 n+1a A: ' Tim ' Ax ' Ax ’ Ax N2Qg_fl: (not nearby the interface) 385: 3833 ' 2TN-1,m+1 sz + TN,m+l(1 + 2 2) ' T Ax Nm (8.22) (8.23) (8.24) (8.25) 138 8.5 EXAMPLE PROBLEM A computer program MELT has been written based on this method to solve the one-dimensional melting problem. This computer program is tested on an example whose analytical solution is available. According to Example 10-2, Ozisik, 1979, the analytical solution of the interface location for a semi-infinite body with constant surface temperature at x - O is u(t) - 2t1/2 (8.26) Although the problem geometry described in Figure 8.1 is a finite slab, both the semi-infinite body and finite slab should behave the same at early times. Again, since the problem in Appendix 8.2 has a heat flux boundary condition an equivalent heat flux can be derived from the constant temperature boundary condition used in Example 10-2, Ozisik, 1979. This heat flux boundary condition then is used to calculate the numerical results. Comparison of the results are shown in Figure 8.2. They display significant agreement between numerical and analytical results. A slightly lagging in the numerical values is observed which is due to the backward difference approximation. The error can be reduced by taking smaller At's. The surface heat flux derived from constant temperature boundary condition goes to infinity at time zero. This might also have contributed to the errors. 139 muaamom aeomuhfiwn< was ~00~HOI52 mo confluealoo N.‘ shaman Amy oEfi omd m ..o N P .0 mod sod 00.0 - n In — P b b h n r n — n n b — n n h o .. Bomboco ..l. * o o .. 32.5.55 0 . .. w: H. - .....IN.O a. .. o r U. . D l O 1 a I e e I a ... . . Ted 1 r o I e e e I m. It . . rl e m. .. rm o u H . . . . ... M7 I D D ‘ D U T. 0 ..l e e e e Im.o ) I! I w w I T mw . T ( r . . _ _ _ u . _ _ . . _ _ t _ P . . t O _. LIST 01' REFERENCES LIST 0? REFERENCES Alifanov, O.M., 1973, Regularization of Solutions of Inverse Problems of Heat Conduction, Heat Transfer - Soviet Research, v 5 n 4, pp. 163-169. Alifanov, O.M. and Artyukhin, E.A., 1975, Regularized Numerical Solution of Nonlinear Inverse Heat-Conduction Problem, J. Engineering Physics, v 29 n 1, pp. 934-938. Alifanov, O.M. and Kerov, N.V., 1981, Determination of External Thermal Load Parameters by Solving the Two-Dimensional Inverse Heat- Conduction Problem, J. Engineering Physics, v 41 n 4, pp. 1049-1053. AlifanOV. O.M.. 1986. MW £1y1ng_ygh1glg§, Mashinestroenie Publishing Agency, Moscow. Alifanov, O.M., Tryanin, A.P., and Lozhkin, A.L., 1987, Experimental Investigation of the Method of Determining the Internal Heat- Transfer Coefficient in a Porous Body from the Solution of the Inverse Problem, J. Engineering Physics, v 52 n 3, pp. 340-346. Bass, 8.8. and Ott, L.J., 1980, Finite Element Formulation of the Two- Dimensional Nonlinear Inverse Heat Conduction Problem, International Computation Technology Conference, ASME Century 2 - Emerging Technology Conference, v 2, San Francisco, California. Bass, B.R., Drake, J.8., Ott, L.J., 1980, ORMDIN: A Finite Element Program for Two-Dimensional Nonlinear Inverse Heat-Conduction Analysis, NUREG/CR-l709, ORNL/NUREG/CSD/TM-l7, U.S. Nuclear Regulatory Commission, Washington, D.C. Beck, J.V. , 1970, Nonlinear Estimation Applied to the Nonlinear Heat 140 141 Conduction Problem, Int. J. Heat Mass Transfer, v 13, pp. 703-716. Beck, J. V. and Arnold, K.J., 1977, Parameter Estimation in Engineering and Science, John Wiley & Sons, Inc., N.Y., N.Y. Beck, J.V., Litkouhi, 8., and St. Clair, C.R. Jr., 1982, Efficient Sequential Solution of the Nonlinear Inverse Heat Conduction Problem, Numerical Heat Transfer, v 5 n 3, pp. 275-286. Beck, J.V. and Murio, D., 1984, Combined Function Specification Regularization Procedure for Solution of Inverse Heat Conduction Problem, AIAA Paper No.AIAA-84-O49l. Beck, J.V. and Blackwell, B. and St. Clair, C. R. Jr., 1985.1, Inverse flgg5_Q2nggg§1g3_111;£2§gg_£;g§lgm, Wiley and Sons, Inc., N.Y., N.Y. Beck, J.V., Blackwell, B.F., and St. Clair, C.R. Jr., 1985.2, Sequential Regularization Solution of the Inverse Heat Conduction Problem, Numerical Methods in Thermal Problems, Proceedings of the Fourth International Conference, Swansea, Wales. Beck, J.V. and Tu, 8., 1988, IHCPZD: Computer Program for Solution of General Two-Dimensional Inverse Heat Conduction Problems, Technical Report for Acurex Corp., Dayton, Ohio, September 1988. 8e11, G.E., 1985, Inverse Formulations as a Method for Solving Certain Melting and Freezing Problems, Numerical Methods in Heat Transfer, Seattle, Washington, John Wiley & Sons, New York, New York. Blackwell, B.F., 1981, An Efficient Technique for the Numerical Solution of the One-Dimensional Inverse Problem of Heat Conduction, Numerical Heat Transfer, v 4, pp. 229-239. Bobula, E. and Twardowska, K., 1985, On a Certain Inverse Stefan Problem, Bull. Pol. Acad. Sci. Tech. Sci., v 33 n 7-8, pp. 359-370. Bonacina, C., Comini, G., Fasano, A., and Primicerio, M., 1973, Numerical Solution of Phase-change Problems, Int. J. Heat Mass Transfer, v 16, pp. 1825-1832. 142 Busby, H.R. and Trujillo O.M., 1985, Numerical Solution to a Two- Dimensional Inverse Heat Conduction Problem, Int. J. for Numerical Methods in Engineering, v 21 n 2, pp. 349-359. Conte, S.D., and de Boor, C., 1980, e ume ' a Analysis - An Algg;1§hm19_5nn;9§gh, McGraw-Hill, Inc. Carnahan, B., Luther, H.A., and Wilkes, J.O., 1969, Appligg_flumggiggl Mgthggs, Wiley & Sons, Inc., N.Y. Dogru, A.H. and Seinfeld, J.H., 1981, Comparison of Sensitivity Coefficient Calculation Methods in Automatic History Matching, Society Petroleum Engineering J., October, pp. 551-557. Flach, G.P. and Ozisik, M.N., 1986, Whole Domain Function Specification Method for Linear Inverse Heat Conduction, Annals of Nuclear Energy, v 13 n 6, pp. 325-336. Frank, I., 1963, An Application of Least Square Methods to the Solution of the Inverse Problem of Heat Conduction, J. Heat Transfer, v 85, pp. 378-379. Goodrich, L.E., 1978, Efficient Numerical Technique for One-Dimensional Thermal Problems with Phase Change, Int. J. Heat Mass Transfer, v 21, pp. 615-621. Haug, E.J. and Arora, J.8., 1979, Applied Optimal Design, John Wiley and Sons, New York, New York. Hensel, E., 1986, Multi-Qimgnsignal Inverse Heat Conduction, Ph.D. Dissertation, New Mexico State University, New Mexico. Hensel, E. and Hills R.G., 1987, Steady-State Two-Dimensional Inverse Heat Conduction, Proceedings, the Sgggng International Conference on v i eerin Sciences, pp. 267-281, Penn State University, PA. Hills, R.G., 1987, ESTIM: A Parameter Estimation Computer Program, Contractor Report, SAN087-7063. 143 Imber, M., 1974, Two-Dimensional Inverse Problem in Heat Conduction, Int. Heat Transfer Conference, 5th, Proc., Tokyo, Japan, v 1, Pap. Cu2. 2, pp. 174-178. Imber, M., 1975, Two-Dimensional Inverse Conduction Problem - Further Observations, AIAA J., v 13 n 1, pp. 114-115. Imber, M., 1979, Nonlinear Heat Transfer in Planar Solids: Direct and Inverse Applications, AIAA J., v 17, pp. 204-212. Katz, M.A. and Rubinsky, B., 1984, Inverse Finite-Element Technique to Determine the Change of Phase Interface Location in One-Dimensional Melting Problems, Numerical Heat Transfer, v 7 n 3, pp. 269-283. Kerov, N.V., 1983, Solution of the Two-Dimensional Inverse Heat Conduction Problem in a Cylindrical Coordinate System, J. Engineering Physics, v 45 n 5, pp. 1245-1249. Lazuchenkov, N.M. and Shmukin, A.A., 1981, Inverse Boundary-Value Problem of Heat Conduction for a Two-Dimensional Domain, J. Engineering Physics, v 40 n 2, pp. 223-228. Lee, T., Kravaris, C., and Seinfeld, J.H., 1986, History Matching by Spline Approximation and Regularization in Single-Phase Areal Reservoirs, SPE Reservoir Engineering, September, pp. 521-534. Markin, A.D. and Pyatyshkin, G.G., 1987, Regularization of the Solution of the Inverse Heat-Conduction Problem in a Variational Formulation, J. Engineering Physics, v 52 n 6, pp. 712-721. Mulholland, G.P., Gupta, P.B., and San Martin, R.L., 1975, Inverse Problem of Heat Conduction in Composite Media, ASME Paper, 75-WA/HT- 83. Osman, A.M. and Beck, J.V., 1987, Nonlinear Inverse Problem for the Estimation of Time-and-Space Dependent Heat Transfer Coefficients, AIAA Paper No. AIAA-87-0150. Ozisik, M.N., 1979, Heat Conduction, John Wiley & Sons, Inc., N.Y., 144 N.Y. Potter, M. C., and Foss, J. F., 1975, Fluid Mechanics, Great Lakes Press, Inc. Romanovskii, H.R., 1984, Solution of Inverse Problems with an Unknown Model of the Process, J. Engineering Physics, v 45 n 3, pp. 1076- 1082. Shapiro, A.8., 1984, TOPAZ - A Finite Element Heat Conduction Code for Analyzing 2-D Solids, Lawrence Livermore National Laboratory, Livermore, CA. Stolz, G. Jr., 1960, Numerical Solution to an Inverse Problem of Heat Conduction for Simple Shape, J. Heat Transfer, v 82, pp. 20-26. Szczurek, J., 1979, Some Approach to the One-Dimensional Stefan Problems (Direct and Inverse), Bulletin de l'Academie Polonaise des Sciences, Serie des Sciences Techniques, v 27 n 2, pp. 215-222. Tang, L.J., 1985, Numerical Algorithm for Inverse Problem of Two- Dimensional Heat Conduction Equation, Numerical Methods in Thermal Problems, Proceedings of the Fourth International Conference, Swansea, Wales. Tikhonov, A. N. and Arsenin, V. Y., 1977, ution o 11- osed Ezghlgmg, V. H. Winston and Sons, Washington, D.C. Tryanin, A.P., 1987, Determination of Heat Transfer Coefficients at the Inlet into a Porous Body and inside it by Solving the Inverse Problem, J. Engineering Physics, v 52 n 3, pp. 346-351. Tu, S. and Beck, J.V., 1987, Sequential Function Specification Approach for One-Dimensional Inverse Melting Problem, Proceedings, Second e o o c o v ce t n En ineerin fing_§21§n§2§, pp. 283, Penn State University, PA. Twomey, S., 1963, On the Numerical Solution of Fredholm Integral Equations of the First Kind by the Inverse of the Linear System 145 Produced by Quadrature, J. Assoc. Comp. Mach. 10, pp. 97-101. Twomey, S., 1965, The Application of Numerical Filtering to the Solution of Integral Equations Encountered in Indirect Sensing Measurements, J. Franklin Inst., v 279, pp. 95-109. Voller, V. and Cross, M., 1981, Accurate Solutions of Moving Boundary Problems Using the Enthalpy Method, Int. J. Heat Mass Transfer, v 24, pp. 545-556. William, S.D. and Curry, D.M., 1977, An Analytical Experimental Study for Surface Heat Flux Determination, J. Spacecraft, v 14, pp. 632- 637.