. ....r. 4 V0 4.... ...,....:..4'.V..a... o... if?! ”WV-d: I!!! 6-! “UL . .ni .J 5 u. E .. .l.-I..x..oz fo‘.). :1! I. .1... an .1.:fl.ur..,.,f_. fit... ...I. p I 4 VA‘Y ,4 o a ’. 'V: .- ‘4 lfr'z I. - ,1:- ‘ r a . ‘I ( . 'l . .3...) .1: l}:-' l I .nl‘tufa.- . '1 .I..m.:r.(l .11. 1P4 J?! lrIavrrffvt...il"IrV-.Jr, v4! a?! f‘ ’ n I I v v . .u....... 4 u. an 3:. r {It '71:! 5.. roil- Alvlv . [ti-H.I.IIVJ.QIV I/Vo': 71? I 9 {Um T I! II‘. ..!..¢ ’43. . a: B. M. ,. N g, .13, r SN 5 , is . CHI E . Q 93 1 mp M m m CR DES . — 4 V . .. .... . .. . V . . .. y. . V , . .7. v...:... II: 2.1:? Ilrlr ”30.145“: I .. l.dr:11~fl.hr.hrt.7.,r||.. . 3., , . . . . .w . K , . . A ‘ .l: I y I ‘ I .A.r . . .l. 7.. . ,( . J»: .V. . . . 1:57.25 S I. ..¢_->.,. . .1 (.3... .I _ t)......| ...__:a.:I-yl tkIBIkL. .Nc . .rf. 35.1.3. .. .th 7 . a. 1 If! 0‘4": I III. v1... 9... . LL... 4.5%.... $5“. 4 1.1.9225... _ . .Va.§§.né...1.§!. . 1 V + L I BRAR Y M ich ' a?” n was)?” ABSTRACT MODEL BUILDING INCORPORATING DISCRIMINATION BETWEEN RIVAL MATHEMATICAL MODELS IN HEAT TRANSFER BY Gerald James Van Fossen, Jr. The primary concern of this thesis is the develOp- ment of mathematical models from experimental data in heat transfer. A model building procedure is prOposed which enables an experimenter to deve10p several mechanis- tic, mathematical models that describe a complex phenomenon. A discrimination criterion is incorporated in the model building procedure to help an experimenter decide which of several rival models describes the phenomenon best. The model building procedure involves performing Optimum experiments to estimate the parameters involved in prOposed models. A criterion for finding optimum experiments for parameter estimation is discussed. Optimum experimental conditions for estimating the parameters involved in a math- ematical model that describes the temperature distribution in a melting solid are found. To illustrate the procedure proposed for model build- ing including discrimination, a mathematical model is found for the specific case of a melting solid. An Optimum, transient experiment for parameter estimation is performed on a finite, one-dimensional, melting, low temperature ufi evaenciituz can be deg Gerald James Van Fossen, Jr. alloy. The temperature data from this experiment is used to find a mathematical model for the melting solid. The proposed model building procedure is generally applicable to experiments for which, a) large amounts of data are generated for each experiment, b) the model is unknown, c) the model can be solved numerically with reasonable expenditures of computer time and d) optimum experiments can be designed and run. WT? (Hy -— h-r ain‘t. I'll ‘ll in MODEL BUILDING INCORPORATING DISCRIMINATION BETWEEN RIVAL MATHEMATICAL MODELS IN HEAT TRANSFER By Gerald James Van Fossen, Jr. A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1973 CHAPTER I DESCRIPTION OF THE PROBLEM. 1.1 1.2 3.79- "I? - '5 is he? H \N H A 1.6 1.7 CHAPTER 2.1 2.2 2.3 2.4 CHAPTER 3.1 3.2 3.3 TABLE OF CONTENTS IntrOduCtionooo00000000000000... 0000000000000 Definition of model building and descrimination between rival mathematical mOdelSOOOOOOOOOOO Possible applications of model building and discrimination between rival models...... Parameter estimation and its use in model building and discriminationOCOOO0.0.0....O... Background of model building and discrimination in other fields.... Application to heat transfer....... The proposed method................ II OPTIMUM EXPERIMENTS FOR PARAMETER ESTIMATION.. Introduction...... Sensitivity coefficients........... The criteria for experiments...... design of Optimum Results for several sample cases............. III EXPERIMENTAL EQUIPMENT AND IECHNIQUE.... The data aquisition system......... Equipment and test specimens.... ...... ....... The experiment.... ii 19 21 28 28 29 31 36 69 69 72 76 R l 2 I) R D. and O o o .Y_~ AL A‘ My. 4 4 4 an... 1A «\w . "I v Mm Db L APPLICIX CHAPTER IV TEST OF THE PROPOSED METHOD............. 4.1 IntrOductionO....0.00000000000000000000000.00 4.2 The model building process................... 4.3 Discrimination between rival mathematical models.. ..... ................................ CHAPTER V CONCLUSIONS .............................. LIST OF REFERENCES ......... ............. ....... .. ...... APPENDIX A FINITE DIFFERENCE METHOD FOR CALCULATING TEMPERATURE DISTRIBUTIONS IN ONE DIMENSIONAL FREEZING-MELTING AND TEMPERATURE VARIABLE PROPERTY PROBLEMS.. A.l General finite difference method............. A.2 The solution of freezing-melting problems.... APPENDIX B PARAMETER ESTIMATION PROCEDURE AND THE CALCULATION OF SENSITIVITY COEFFICIENTS. B.l The Gauss-Newton or linearization method of minimizing the sum of squares function.... B.2 Calculation of sensitivity coefficients...... APPENDIX C THE CALCULATION OF SURFACE HEAT FLUX FROM AN INTERNAL TEMPERATURE HISTORY.... iii 81 81 81 102 106 109 112 112 119 135 135 138 142 TABLE 2.4.1 2.4.2 2.4.3 4.3.1 LIST OF TABLES + KOpt for several values of TT with pL the parameter to be estimated for the semi-infinite exampleooooooooooooooooo 25‘ + + opt for several values of TT and LT with T and pL being the parameters to be estima? Locations and times for maxima and minima of sensitivity coefficients figL and 8% ... f Comparison of the likelihood functions of the four models...................... iv ed for the semi-infinite example.... 44 49 67 104 1' ‘[[[[[([[‘l‘||‘lllll I a .0 :w ( /l\ 10. l 1 an. .... 2 q. 1L 1 2 «In FIGURE 1.2.1 1.7.1 2.4.l(a) 2.4.l(b) 2.4.2 2.4.3 2.4.4 2.4.5 2.4.6 LIST OF FIGURES PAGE Concentration of B versus reaction time for the two models A-DB“...C and A-O-B-O-C.... Flow diagram of the prOposed model building procedure............................ Dimensionless temperature of a freezing semi-infinite body subjected to a step decrease in temperature at the boundary for several values of the dimensionless heat of fusion, L$.................... ..... ... Dimensionless temperature of a freezing semi—infinite body subjected to a step decrease in temperature at the boundary for several values of dimensionless fusion temperature, T;................. ...... . Dimensionless sensitivity coefficients for a freezing liquid semi-infinite body with a step decrease in temperature at the surface ............ ............ ............... [Sept versus dimensionless heat of fusion for a single thermocouple with PL the property to be estimated............ ...... .... Merit function, M, versus £5Tmax"'°""°’°"° Product of merit function, M, and ZSOpt versus [5T used to determine the max optimum AT for estimating PL........ ..... max Sensitivity coefficients at x/E = 0.0 for a melting-freezing,finite body with constant heat flux, q, at x=O and insulated at x=E ..... 23 40 41 42 45 45 47 50 FIGURE 2.4.7 2.4.8 2.4.9 2.4.10 2.4.11 2.4.12 2.4.13 2.4.14 2.4.15 2.4.16 PAGE Sensitivity coefficients at x/E = 0.5 for a melting-freezing, finite body with constant heat flux, q, at x=O and insulated at =EOOOOOOOOOOOOOOOO......OOOOOOOOOOOIOO ..... .0 Sensitivity coefficients at x/E = 1.0 for a melting-freezing, finite body with constant heat flux, q, at x=O and insulated at x=E................ .............. Dimensionless, 33L, versus dimensionless time,‘T , for a finite body with constant heat flux, q, at x/E= O and insulated at x H=l.........OOOOQOOOOOOOOOOO0.0......0.... E versus Tm for several values of ATmax for a melting, finite body of solid initially at the fusion temperature with PL ' the parameter to be estimated................. Egopt versus lkTmax for finite body of solid initially at the fusion temperature with pL the parameter to be estimated................. The product Mrsop t versus csTmax used to determine the optimum [5T for estimating max the para-mater PLOOOOOOOOOOOOOOOOOOOOOOOOOOOO Experimental heat flux from c0pper heater- calorimeter and heat flux through copper lid versus time............................... PL and STf at x: O. 0 versus time for a finite body insulated at x=E and the heat flux of Figure 2. 4.13 prescribed at x: O. O....... ...... . . .. ..... ... Sensitivity coefficients ST Sensitivity coefficients SabL and S x/E = 0.5 versus time for a finite body insulated at x=F and the heat flux of Figure 2.4.13 prescribed at x=0.0 ............ . Sensitivity coefficients 33L and s; f x/E = 1.0 versus time for a finite boay insulated at x=E and the heat flux of Figure 2.4.13 prescribed at x=0.0. .......... .. at vi 51 52 55 58 59 61 63 64 65 66 FIGURE 3.2.1 3.3.1 4.2.1 4.2.2 4.2.3 4.2.11 4.2.12 4.2.13 The experimental equipment............. ......... 74 Partial section view of nylon cup specimen showing copper lid in place and location of thermocoupleSOOOOOOOOO0.0.0.0.........OOOOOOOOOO77 Constant thermal conductivity and heat capacity model and the data versus time.................. 83 Residuals for constant thermal conductivity and heat capacity model versus time.. ....... .... 84 Residuals for the constant thermal conductivity and heat capacity model versus measured temperature.....OOOCOOOOOOOOOOOOOO0.00.......... 85 Step change in thermal properties model and data versus time ..... ....... ..... .......... ..... 87 Residuals for the step change in thermal property model versus time...................... 89 Residuals for the step change in thermal properties model versus measured temperature.... 90 Energy absorbed by melting over a range of temperature model........................ ....... 92 Melting over a range of temperature model and the data versus time...................... ...... 93 Residuals for the melting over a range of temperature model versus time................... 95 Residuals for the melting over a range of temperature model versus measured temperature... 97 Isothermal melting-freezing model and the data versus timeOOOOOOOOOOOOOO0.0.0....00.0.0... 98 Residuals for the isothermal melting-freezing model versus time....................... ........ 99 Residuals for the isothermal melting-freezing model versus measured temperature.. ............. 100 vii mm A.) 5.1 c 2 A.l A 2 M A) ‘- I‘ ‘II ‘ ‘ . [3‘ ( [ [ ( I l [ I ll ‘ l ‘7‘ (2.2 FIGUR A.1 A.2 A.3 A.4 A.5 B.1 C.1 0.2 0.3 PAGE Node geometry and nomenclature for the finite difference solution to the heat conduction equation at the interface between dissimilar materials.......... ...... ....114 Node geometry and nomenclature for the finite difference solution to the heat conduction equation at the heated surface with a known heat flux..... ..... .............................114 Node geometry and nomenclature for the floating node used in the finite difference solution of the freezing-melting problem........114 A comparison of the finite difference solution to an exact solution of the freezing-melting prOblemoooooooooooo00000ooooooooooooooooooooooo0131 A comparison of the finite difference solution to an approximate solution of the freezing- melting prOblemOOOOOOOOOOOOOOOOO0.000.000.000000134 A comparison of exact and numerical approxima- tion of the dimensionless sensitivity co- efficient Si (x,t) for a finite body with constant heat flux, q, at x = O and insulated atx:EOOOOOOOOOOOOOI.......OOOOOOOOOO...0.0.0.140 A comparison of exact and numerical approxima- tion of the dimensionless sensitivity co- efficients sT (x,t) for a finite body with constant heat lux, q, at x = O and insulated atx:EOOOIIOOOO.........-OI0.00.00.00.000000000141 Ratio of temperature rise to maximum temperature rise versus ratio of dimensionless time to dimensionless duration of heat pulse at the insulated surface of a plate with heat applied to the other side. ..... ................. . ..... 144 Illustration of how a time varying heat flux is broken up into a series of heat pulses. ..... .144 Comparison of heat flux calculated from temperature data accurate to six places for 1, 2, and 3 future temperatures for a slab with a step in flux at one surface and the Other surface insulated........... .............. 147 viii FIGURE PAGE C.4 Heat flux calculated from noisy data for 2 future temperatures for a slab with a step in flux at one surface and insulated at the other surface................. ..... .........148 ix 03> (n). D(n) CNBD erf(x) erfc(x) L LIST OF SYMBOLS (n), B(n), Coefficients used in numerical approximation to heat con ucti n equation defined by equations A.1.% Finite difference parameter defined in Appendix A Specific heat at constant pressure Length of finite body x 2 .2. -t Error function erf(x) =-Vfif o 6 dt Complementary error function erfc(x) = l-erf(x) Constant heat flux applied to surface of a semi-infinite body Thermal conductivity Ratio of thermal conductivity of liquid phase to thermal conductivity of solid phase Heat of fusion L(model i/data) Likelihood function for model 1 given + + LT’ Lq M m N n P 9 the data Heat of fusion non-dimensionalized with respect to temperature and heat flux respectively Merit function for maximum temperature rise Number of discreet times in data Number of data points number of thermocouples Vector of parameters with elements pi Heat flux Mm) |N IR: Sensitivity coefficient — derivative of temperature T with respect to parameter p1 Initial temperature Fusion temperature Surface temperature (x = 0) Fusion temperature non-dimensionalized with respect to temperature and heat flux respectively Dimensionless temperature of solid and liquid phase respectively Temperature at location 1 and time t Time Vector of predictions with elements Wu Location Vector of observations with elements yfi Sensitivity matrix with elements 213 xi €\A£.PUTV{ 9+9 0 ’d d- D Dll>| [DH 53 N I'Tb >’ I“ macro GREEK SYMBOLS Thermal diffusivity Ratio of thermal diffusivity of solid phase to thermal diffusivity of liquid phase Matrix formed by the product g'; Matrix whose elements are égj defined by equation [2.3.2] Determinant of the matrix 4; Optimum Kas defined by equation [2.3.6] Maximum temperature rise experienced during an experiment Location of solid-liquid interface Root of equation [2.4.4] or finite difference parameter defined n Appendix A Vector of controllable variables Density Variance of observations Dimensionless time Weighting matrix used in likelihood function xii CHAPTER I DESCRIPTION OF THE PROBLEM 1.1 Introduction One of the most fundamental problems in science is that of obtaining mathematical models from observations. The basic objective of this thesis is to propose and demonstrate a procedure that can be used to develop an adequate mathematical model for certain types of physical PI‘Ocesses. The procedure is successfully applied to the Specific case of finding a mathematical model from experi- In‘E‘I‘I‘t:al data for a melting, one-dimensional solid which is heated at one surface and insulated at the other. This re- presents a complex heat transfer process. The proposed pr"Deedure could be used to find an adequate mathematical InOC1621 for any process where one can obtain large amounts of data from a single transient experiment at a relatively low cost . The individual parts of the prOposed procedure are not new; however, the idea of using optimum experiments and discrimination coupled with actual heat transfer data for find ing a mathematical model is original. The prOposed process is called model building which incOrporates discrimination between rival mathematical m . . 0C‘eels, these terms are defined in section 1.2 of this ‘ 1 I I chapter which c some ex in heat paranet ature c emetica; describe discrim; .2 chapter. Model building procedures have been suggested which do not include discrimination. Section 1.3 contains some examples of how the prOposed procedure would be of use in heat transfer. Section 1.4 discusses the importance of parameter estimation in model building. A review of liter- ature concerned with model building between rival math- ematical models is contained in section 1.5. Section 1.6 describes how the existing methods for model building and discrimination can be applied to heat transfer. The pro- posed procedure for developing an adequate mathematical. model is given in section 1.? along with a flow chart that outlines the proposed procedure. 1. 2 Definition of model building and discrimination between rival mathematical models. Model Building Model building is the complete procedure which en- ables an experimenter to develop a "best" model for the pro cess. It is not simply the derivation of mathematical models from conservation principles. All the models in this thesis can be so derived. The experimenter may; however, be able to derive several of these models without knowing Which is the correct one. In the model building procedure the simplest reason- 8L"ale "first cut" model is proposed and fitted to the experi- mental data (i.e., the parameters in the model are estimated). 1‘ «e quality of fit is then examined to learn the nature of the deficiencies in the model. based on this knowledge, 3 a modified model is proposed to better describe the process. Discrimination A very important part of model building is called discrimination. Discrimination involves the design of ex- periments so that differences between the predicted responses of two or more models is the greatest. By using such ex- periments along with a decision-making criteria, one can better decide which is the correct model. Various criteria are presented in section 1.5. An example of the need for design of experiments for discrimination has been given by chemical engineers (l). A substance "A" reacts in the presence of a catalyst to form substance "B" which in turn forms "C". One possible model for this is given by A-O-B‘C [l . 2 . 1] Another possible model is 4‘320 I [i . 2 . 2] The concentration of "B" versus reaction time is shown in Figure 1.2.1. If the reaction is observed only until time tl, no discrimination can be accomplished because the pre- dicted responses and the data are essentially the same until t1. All the parameters in the model represented by equation [1.2.2] cannot be determined accurately if the reaction is observed only until time t1, because there has been no Significant production of "B" from "C" for Ofitétl. If the J‘i b. A-a-B:=:C -'-'- A-iiB-Iflz ‘7 Data Concentration of B ¢—-————-—-—-—- Time Figure 1.2.1 Concentration of B versus reaction time for the two models A-a—B::=C and A-aiB-ufl}. applied number c One app radiant mid. if aPPlica‘ They US: the Cha; U893 a “0591 b any Sit Reeds I. 5 reaction is allowed to proceed until time t2, we see that the model of equation [1.2.2] is the better of the two proposed. 1.3 Possible application of model building. Model building, including discrimination, can be applied to many areas in heat transfer. For example, a number of models have been proposed for ablating solids. One application to this area is to determine if significant radiant heat transfer exists within the ablating material and, if present, should be included in the model. Another application relates to the work of Pfahl and Litchel (2). They used an isothermal change of phase model to describe‘ the charring of cork but a competing and untested model uses a chemical reaction to describe the charring process. Model building, including discrimination, can be used in any situation where the mathematical model is unknown or needs to be improved. 1.4 Parameter estimation and its use in model building. Parameter estimation is a discipline that provides techniques for calculating many of the parameters in a math- ematical model from a single complex experiment. Some engineers use the terms nonlinear estimation, identification, nonlinear regression and nonlinear least squares instead of parameter estimation. W47 4' "31" LII" appear: to desi by a s: "ccmple perimen meters. 6 The classical method of determining parameters appearing in a relatively complex mathematical model is to design several experiments, each of which is described by a simple model involving one parameter. Thus if the "complex" model contains two parameters, two separate ex- periments would be conducted in order to find both para- meters. An example is linear heat conduction involving two parameters, thermal conductivity and heat capacity. If a steady state experiment is performed, only the thermal conductivity is involved in the reduced math- ematical model. If a transient experiment is performed for a specimen with no temperature gradients, only the heat capacity is involved in another reduced model. In most cases such separate experiments would be adequate, but there are materials which dry or otherwise change during the prolonged tests uSually involved. When these changes occur, a rapid transient experiment might be possible that would permit the estimation (using parameter estimation) of the properties. snauld ‘ I qua; 143:; other v 1.4.1 Distinction between_parameters and properties The distinction between parameters and properties should be made clear. A property is a characteristic quality of a system and may be a function of any of the other variables in the problem. A parameter need not have the physical significance of a property and it is not a function of any other variable. A parameter may be con- tained in the expression for a property; for example,in the case of temperature dependent thermal conductivity one may write k(T) 2 k3 + kT [1.4.1] where k(T) is the prOperty thermal conductivity and k8 and k are parameters. 1.4.2 Parameter estimation criterion In testing a model's ability to describe a physical phenomenon, the predicted response is compared Ivith experimental data taken from the actual physical saystem. In making this comparison, the parameters in the nusdel can be estimated from the data if they are unknown. A basic tool of parameter estimation is the minimi- za13ion of a weighted sum of squared differences between the measurements and predicted values from the mathematical 8 model. (See Appendix B for a more complete description.) The parameters in the model are varied until the combination is found that minimizes the weighted sum of squares. The ability of the digital computer to solve complex mathemati- cal models has allowed parameter estimation to become more widely used. Box and co-workers (3’4) have develOped the method from a statistical viewpoint. Beck(5’6) was the first to apply the method to finding parameters involved in a partial differential equation in heat transfer. He used the method to find thermal conductivity and heat capacity simultaneously from a single transient experiment. Pfah1(7’8) used the method to find seven parameters at once in charring cork. 1.4.3 Importance of optimum experiments in parameter estima- tion. It is important that parameter estimation experiments be carefully planned in advance. The best or optimum ones should be run in order to obtain the greatest accuracy. One might think that experimentalists always do this. This is not usually true in the sense meant here, however. One can run the best experiment in terms of careful specimen prepara- tion, placement of sensors, test procedure, and data aquisi- tion and still have a poor experiment for simultaneously estimating several parameters. It may be clear how to de- sign an experiment for estimating a single parameter. It is not obvious; however, how to design an Optimum experi- ment when several parameters are to be estimated. Optimum parameter estimation experiments are discussed further in Chapter II. l.§iBackground of model building and discrimination in other fields. Model building is the name of the general procedure that we are discussing. It has a number of components in- cluding selecting optimum experimental designs, performing the experiments, analyzing the data, etc. These are dis- cussed as a whole in Section 1.7. An essential step in model building is discrimination. It is so important a step that there is a tendency to confuse discrimination with the the entire process. Discrimination is only a part of model building, but it is a crucial part in the proposed procedure. 1.5.1 Model building There is extensive literature on the mathematical modeling of different phenomena. There are cases; however, for which the usual continuum mechanics approach can produce several competing models. Model building is needed to choose the "best" one of these rival models. The model building literature is much smaller than the mathematical modeling literature. Box and Hunter Most of the model building work was done by chemical engineers and chemists in cooperation with statisticians. (9) Box and Hunter proposed a method for use in model build- ing. They suggested treating the estimated parameters in the 10 model as observations and the controllable variables as in- dependent variables. Controllable variables are the varia- bles an experimenter is able to change or control. If the model is correct, the estimated parameters should not change as the settings of the controllable variables change, i.e., parameters should stay constant. Factorial design1 is used to choose the setting of the controllable variables (the fac- tors) and the parameters are estimated for each setting of the factors. Statistical analysis is used to determine the effect of the factors on the parameters thus giving the ex- perimenter important clues as to how the model should be mod- ified. In another paper Box and Hunter(1) state "all infor- mation relating to the possible inadequacy of a tentatively entertained model is contained in the residuals". Residuals are the difference between the measured response and the re- sponse predicted by the model. Examination of plots of the residuals can aid one in telling where improvements should be made. Residuals should be plotted in any way that might shed light on pertinent questions. Anscombesll) Anscombe and Tukeyslz) and Draper and SmithSlB) in addition to Box and Hunter, all recommend graphical examination of the residuals to show deficiencies in the model. 1.5.2 Discrimination Discrimination has several aspects. These include, .__‘ lFactorial design is a statistical procedure that allows one to determine efficiently the effect of the factors on the response; see reference (10) page 122. 11 a) the develOpment of a discrimination criterion, b) the use of this criterion and c) the optimum design of experiments to promote the discrimination. Parts a) and c) can be confused because the optimum experiment design can involve a criterion based upon the discrimination criterion. Discrimination between rival mathematical models is an area known to statisticians as hypothesis testing. Work has been done by statisticians and engineers on both the discrimination criterion and on designing the most efficient experiment to accomplish the discrimination. Hunter and Reiner Hunter and Reiner(14) proposed a design scheme for the case of two rival models. Their experimental design consists of making n initial runs to estimate the parameters in each model (n must be greater than the number of para- meters in each model). The n+1St setting of the controlla- ble variables is chosen so that the predicted responses of the two models are farthest apart. The n+18t data point is taken, the parameters estimated again, and the n+2nd setting of the controllable variables determined in the same manner. The procedure is thus sequential. The discrimination is accomplished by making replicate runs to get an estimate of the variance which is a measure of the experimental error. The sum of squares of the residuals compared to the variance in an F-test is the discrimination criterion. An F-test is a statistical test used to guard against making the wrong decision by comparing the residual sum of squares with an estimate of experimental error. mu? 0. O YCEI tor 1985 I QM +5 ..62‘ 4016 ECCC -‘. 3x.» A» ...» at. is e 2» 5v ‘1‘ . s in .MM .3 r .» ..\ .v «\V Ali. ‘1‘ “Us \‘ IRE .“ |‘\ 12 Roth Roth(15) extended the experiment design criterion of Hunter and Reiner to any number of rival models and used a Bayesian approach for the discrimination. Roth defined 0 = lwgl)-w§i)uw§2)-w§i)uw§3)-wgi) 13 . ... [1.5.1] where wgi) is the predicted response of model i with the vec- tor of controllable variables set at Ej‘ Cij is a direct measure of spread between the predicted responses. Roth chose to select‘gj, the setting of the controllable variables, to subject the current "best" model to the severest test; to accomplish this he defined the utility function Utgj) = EE:P(model 1)ocij [1.5.2] 1 where P(model i) is the probability or degree of belief in model i. The experiment desigg procedure is to select the ‘53 that maximizes Ugéj). Roth used Bayes' theorem to combine a prior probabil- ity for each model and the likelihood function for each model to find a posterior probability for each model; i.e., P(model i/data) = P(model i)-L(model i/data)/P(data)[l.5.3] The maximization ofthua posterior probability is Roth's discrimination criterion. P(model i/data) is the posterior probability of model i or simply the degree of belief in model i in light of the data; P(model i) is the prior proba- bility for model i; L(model i/data) is the likelihood funct- ion for model i which will be discussed shortly; and P(data) is given by 13 P(data) == §E1P(mode1 i)-L(model i/data) [1.5.4] 1 If the errors in the measurements are assumed in- dependent and normally distributed with zero mean and con- stant variance,CTg, the likelihood function can be expressed as N N L(model i/data) = 1 exp - 1 ii: (y -w )2 «is 2? u u “=1 [1.5.5] Where N is the number of measurements, yu are the measure- ments and Wu are the predicted responses. If no prior information is available P(model i) can be set to l/m for all i where m is the number of rival models. The decision of calling one model best is made by comparing the magnitudes of the posterior probabilities for different models. The posterior probability not only Shelps choose between models but gives the probability of eeach model in light of the data. This is also a sequential procedure . Box and Hill The experiment design procedures of Hunter and Resiner or Roth do not take into account the magnitude of tile experimental error; thus, it could be possible for the neaxt design point to be in a region where the confidence regions for the predicted response for each model overlapped a j (l- jill " I l4 and no real discrimination could be accomplished. A con- fidence region is a band around the predicted response from a model where one has a certain degree of confidence that, if the model is correct, measurements will fall within this band. (16) used the concept of entrOpy from Box and Hill communication theory to develop an experiment design pro- cedure to correct this fault. EntrOpy can be defined as m s = - 275an: [1.5.6] i=1 where m is the number of rival models and‘fl; is the prob- ability associated with model i. The least possible information regarding discrimination occurs when entrOpy is maximum; it can be shown that this happens when 1ri = l/m [1.5.7] for all i. The criterion of Box and Hill is‘then to znaximize the expected change in entropy between input and (output; that is, to perform the experiment which will cause 'the greatest change in entropy. An upper bound for the eaxpected change in entropy was given by Box and Hill as m m (C72+(72) .122. . .. i'n'l 3'“’l cr2+cri)«r2+cr3) i=1 j=i+l J + (yn(i)-y£j))2 1 + 1 [1,5,8] (62+of)(oj+o§) 15 "i , n-1 and Trj , n-l i and 3 respectively after n-l runs have been performed, are the prior probabilities of models 0'2 is the constant variance of the data, 0'1 and 05 are (i)‘ the variances under models 1 and j respectively, and yn and y(j) are the predicted responses from models 1 and j n th or next data point. This ex- (i)_ n to the criteria of both Hunter and Reiner and Roth; that is, respectively for the n pression contains the term (y yé'fl)2 which is similar the next data point is taken at the setting of the control- lable variables where the predicted responses are the farthest apart. The Box and Hill criterion also includes the variance of the data and of each model. Hunter and Hill propose exactly the same criterion but they use repeated runs to obtain an estimate of the variance,<72. The discrimination criterion used by Box and Hill is the same as that of Roth, i.e., the Bayesian posterior probability for each model. Meeter, Pirie, and Elot(18) compared the method of Box and Hill with a somewhat similar method developed by Chernoff(19). They favor the Chernoff procedure because the Box and Hill criterion requires maximization of an upper bound; it would seem more logical to maximize a lower bound but this is not easily obtained. Results on several test cases done by heeter, Pirie and Blot; however, show the box and Hill procedure superior. Chernoff's criterion was develOped assuming the cost of making measurements was zero. This is not the case in most engineering problems Where the cost of experiments is generally high. 16 Hunter, Hill and Wichern Hunter, Hill and Wichern<2o) presented a joint design criterion for model discrimination and parameter es- timation. Their experiment design criterion is to maximize C = wlD + w2E [1.5.9] with respect to the controllable variables. Here wl and w2 are non-negative weights, D is essentially the Box and Hill criterion described above, and B is a parameter estimation criterion. Initially wlD is made larger but as experiment- ation progresses and one model begins to emerge superior w2E begins to dominate. The discrimination criterion is the Bayesian posterior probability for each model. Atkinson (21) gave a discrimination criterion which Atkinson is a statistical test for determining if there is any evi- dence that the models give significantly different fits to the data; in other words, it tests the hypothesis that all the models are the same. The test statistic is L? ~ 2 2 HF) (;?§§L [1.5.10] ”(3 I F‘ where N is the number of measurements, p is the number of A models, (YF)2 is defined by A 2 A 2 (m = 2 (5,3,4) [l.5.lO(a)] 17 f‘is the model formed by regression using the observations or measurements as the dependent variable and the predic- tions from each model as the independent variables. ‘? combines the predictions from all the models into a single set of predictions that describe the data. (f§)2 is defined by N ~ 2 __ _~ 2 (m — Z (yj 1’) [151000)] =1 with 3" defined by N Z (fij- ?)2 =1, i=1,2,...,p [l.5.lO(c)] f is the predicted value from model i for data point 3. 13 X, is any constant. Thus 2" is equally distant from each model. (Y?)2 is used as an estimate of error. The test statistic is tested by an F-test with N-p and p-l degrees of freedom to determine if there is any significant difference in the way the models fit the data. Atkinson does not consider.the problem of designing experiments to make discrimination more efficient; his method merely treats each model as a formula that is supposed to represent some aspect of the data and tests to see if they are different. 18 Reilly The ratio of likelihood functions was used by Reilly(22) as his discrimination criterion. The likelihood function can be defined in terms of the probability den- sity function<23). Let f(X,§) be the joint probability density function of the independent, random variables X and parameters 3. Suppose that N observations are made on 1, and let (Y1,Y2,Y3,...,YN) be the observations; then the function given by L(Y1’Y2,Y3,...,YN:§) = f(Yl,g)f(YZ,§)f(Y3,g)...£(rN,§) [1.5.11] is called the likelihood function. A maximum likelihood estimator §,is an estimator that maximizes the likelihood function with respect to the parameters 2. If the errors are independentznui gaussian with zero mean and constant variance, maximum likelihood parameter estimates are the same as those given by least squares. For the purpose of comparing the relative plaus- ibilities of two parameter vectors £1 and 22 it is conven— ient_to compare their likelihoods by examining the ratio. The likelihood ratio L(§flML(§é) is sometimes called the odds ratio and is a direct measure of the relative plaus- ibilities of 3 and 2 in light of the data. According to 1 2 Reilly, a likelihood ratio of 10 is ordinarily taken as showing a real difference while a ratio of 100 shows strong preference for one parameter over the other(22). 19 To discriminate between two rival mathematical models, Reilly suggests that we examine the ratio Ll(§l)/L2(§2) where Ll(§l) and L2(£2) are maximum likelihood functions for models 1 and 2 reSpectively. This gives a comparison of models at their individual best. 1.6 Application to heat transfer. All the methods described above with the exception of the last two were developed for a different type of ex- perimental technique than that used in our heat transfer experiments. In chemical engineering an experiment to study a reaction rate mechanism, for example, may consist of sev- eral runs. A typical experimental run might be to start with a known initial concentration of reactants; fix the temperature, pressure, catlyst, etc. and let the reaction proceed for a certain length of time. The reaction is then stopped and the reactants analysed to find the concentrations of reactant and product. Thus each run may contain only a single data point, and several runs must usually be made before one may estimate all the parameters in the proposed model. In transient heat transfer, however, it may be possible to obtain all the necessary data from a single run. The experiment might consist of heating a sample of the "unknown" material and measuring the temperature as a function of time for several locations on the sample. The temperature measurements are usually obtained with thermo- couples. Several thermocouples may be used which allows 20 several locations at one time to be recorded. The signal from each thermocouple may be recorded continuously or it may be discretized at this point. A number of parameters in the model are estimated from a single experiment. The method for model building proposed by Box and Hunter(l); i.e., treating the estimated parameters in the model as observations, may be misleading in many cases in heat transfer. It is demonstrated in Chapter II that if the settings of the controllable variables are not chosen carefully, the estimated parameters may be inaccurate or it may not be possible to find them all simultaneously. There- fore, this method is not recommended for model building in heat transfer. Many researchers recommend graphical examination of (1’ 11’ 12’ 13) This seems the residuals for model building. to be a most useful tool for model building. By plotting the residuals against different variables one may learn where deficiencies in the model exist and prOpose a new or modified model that will correct them. This process will be illustrated by example in Chapter IV. In situations like the example of the reaction rate study discussed in the beginning of this section, many ex- periments are required to obtain the data. If one wishes to discriminate between two or more models that are supposed to describe the reaction mechanism, some method of conserving (minimizing the number) experiments is needed. For this I?eason the design strategies outlined in section 1.5 were (ieveloped. All the experimental designs developed thus far 21 are concerned with finding the setting of the controllable variables at which to take the Egg: data point for best discrimination, thus all the designs are sequential. In many cases in heat transfer there are a relatively limited number of controllable variables. Frequently all the data necessary for discrimination can be obtained from a single experiment; in such cases sequential designs are not needed. ThefiBayesian approach used by Roth and also Box and Hill allows prior knowledge or belief about each model to enter the design criterion. The experimenter usually has no prior preference for any one model initially, but as experie mentation and sequential discrimination proceed concurrently he may gain evidence that one or more models is superior to the others. This evidence or belief is allowed to enter the decision criterion. Without sequential experimentation, one has no way to gain sufficient a priori knowledge to allow it to enter the discrimination decision. While the Bayesian approach is a very powerful tool, it will not be employed here because of the difficulty of performing different se- quential experiments. On the other hand, both the method for discrimina- tion proposed by Atkinson and the method of maximum likeli- hood ratios described by Reilly seem to be directly appli- cable here. Both these methods involve the residuals in some manner. The method of Reilly is to be employed. 1.7 The proposed model building method. The proposed model building procedure is illustrated 22 in the flow diagram of Figure 1.7.1. If the experimenter lacks sufficient knowledge about the process to propose a complex model, it is recommended that a simple mechanistic model be proposed for the process even though the model is suspected to be deficient. This is block 1. of Figure 1.7.1. in block 2. the Optimum experiment for parameter estimation is found. Estimating the parameters utilizing data from the optimum experiment will cause the model to describe the process as adequately as it can. The optimum experi- ment is then performed (block 3.) and the parameters esti- mated from the data (block 4.). In block 5. the model and the data plotted together are examined visually to try and extract information about the model. The residuals should be plotted against any pertinent variables and examined visually. Valuable clues as to the nature of the defi- ciencies in the model can be gained in this manner. This is a very important step. In block 6. the experimenter must decide if the model or models he has to this point in the process are plausible or capable of describing the features of the process he is interested in. There is no formal criterion for this decision; it is left to the judge- ment of the experimenter. If they are plausible, he pro- ceeds to block 8. the discrimination stage; if they are not, he procedes to block 7. In block 7. the experimenter must evaluate the information he gained in block 5. and propose a new or modified model based on this information. The experimenter then proceeds back to block 2. and repeats the process until he has one or more plausible models. 23 1. First Cut or Guessed Model 1 2. Find Optimum Experiment for _ Parameter Estimation l 3. Perform Experi- ment 7. PrOpose New Model from 4 Fit Model to Information Data (Estimate ggigigifig Parameters) Resid a 5. Analyze Model f and Residuals to Determine Deficiencies 1 6. Are one or more I J Models Judged to No be Plausible? * Yes 8. Perform Discrimination Figure 1.7.1 Flow diagram of the prOposed model building procedure. 24 This procedure will be illustrated in Chapter IV. In many cases the experimenter has some prior know- ledge about the process. If enough information is available to propose one or more mathematical models, then the build- ing of competing models may be skipped entirely and one may go on to discrimination. It is highly recommended that the residuals be examined in any case; they contain all the information about the "goodness of fit" of the models. The next step after proposing new models is to test them to determine which fits the data adequately. This is the discrimination process (block 8. of Figure 1.7.1). Three possible cases arise here: 1.) none of the models is adequate, 2.) one of the models is adequate or 3.) more than one of the models is adequate. Should Case 1.) occur, more models should be developed. Case 2.) is the most desirable case; if this should occur, the investigation can be con- cluded. For Case 3.) two explanations are possible; either one model is a special case of the other or the data is not sufficiently accurate to allow discrimination. If the for- mer is true, the experimenter should be the judge of which model suits his purpose best based on other considerations. If the latter is true, methods of improving the experiment should be sought. To improve discrimination the experimenter might want to devise some special experiment to discriminate between the models. As an example of this, suppose two rival models that describe the temperature history of a solid being heated are : l.) a model that has change of 25 phase and 2.) a model that has a chemical reaction. Since both models might describe the temperature equally well during heating, the experimenter should devise another ex- periment to discriminate. Change of phase is reversible; that is, the specimen can melt during heating and solidify during cooling. Thus for change of phase, the measured temperatures should exhibit the same behavior while heating or cooling. The chemical reaction is irreversible; that is, the reaction occurs during heating but the reactants do not return to their original state during cooling. Hence the measured temperature for the chemical reaction should exhibit considerably different behavior during cooling than it did during heating. With this information the experimenter could discriminate more effectively by performing an experiment that involved both heating and cooling. All statistical tests are based on comparing the re- siduals with some measure of experimental error to determine if the model is adequate in light of the data. The method proposed by Atkinson does not require the experimenter to supply an estimate of the experimental error. This method, however, tests the hypothesis that all the models fit the data equally well; one could just as well say it tests the Ihypothesis that all the models fit the data equally poorly. It provides no means to test how well the models fit the data on an individual basis. The likelihood ratio test requires the experimenter ‘to supply an estimate of the eXperimental error which can sometimes be difficult to obtain. If a reasonable 26 estimate of experimental error is available, this method can be used effectively in discrimination since only the ratio of likelihoods is of concern. The experimenter who is familiar with his equipment can usually supply an esti- mate of experimental error that is good enough for this purpose. The alternative is to perform replicate runs to obtain an estimate of variance. Although replication is generally to be recommended, it will not be used in this work. The likelihood function for zero mean and gaussian error is L(mode1 i/data) = 1,52 EXP [-HX'EYY-J-(X‘Eg (210 [det1!]ig where I is the covariance matrix of the errors. For in- dependent errors 3! reduces to a diagonal matrix with the diagonal elements equal to the square of the variance. If the variance is constant then equation [1.7.1] reduces to [1.5.5] . For correlated errors ‘2: can be a full N x N .matrix. Unfortunately the errors are correlated in a tran- sient experiment involving temperature measurements when the sampling rate is relatively high. This is the type of ex- ‘periment that we have. We have no way to evaluate the elements of the co— variance matrix; all we can do is provide an estimate of the (nonstant variance (experimental error). When errors are correlated it is known that (24.25) 27 N Z (ya—Wu)?» (i-m'yjlq-g) [1.7.2] 1 u: .1. C,-2 Considering this we will define a new likelihood function in an attempt to approximate the true value of the likeli- hood function from [1.7.1] . Let the new likelihood function N . l 2 L(model 1/data) = EXP - E (y -w ) [: 2CT2N u u ] u: [1.7.5] This is not the correct form of the likelihood function but be it is the best we can do at this time. Equation [1.7.3] will be used in the likelihood ratio test to discriminate between rival models. CHAPTER II OPTIMUM EXPERIMENTS FOR PARAMETER ESTIMATION 2.1 Introduction An optimum experiment is defined as the experiment which allows the parameters in a mathematical model to be calculated with the greatest accuracy. The combination of experimental conditions yielding predicted responses that are most sensitive to changes in the parameters are usually the optimim conditions. This will be demonstrated later. If experiments are not carefully planned, the experimental points may be so situated in the space of controllable variables that the calculated parameters may not only be im- precise but also highly correlated(l). Controllable vari— ables are defined to be the variables the eXperimenter has the ability to change; some examples. might be location of sensors, maximum run time for transient experiments, boundary conditions, and initial conditions. If parameters are correlated for certain settings of the controllable variables, they cannot be determined simultaneously. Seinfeld(26) calls this condition the observability problem while others call it the identifiability condition. An example of an experiment with correlated para- meters is a semi-infinite body heated by a constant heat flux 28 29 F0, at the surface x = 0. Here the controllable variable is the location, x, of the measurement. At the heated sur- face the temperature at time t is given by<27) T(O,t) = 2Fo(t/‘fl'p cpk)’1r [2.1.1] where pop is the volumetric heat capacity and k is the thermal conductivity. If a single thermocouple were located at the surface x = O, k and pop can not be found simultane- ously; they are said to be correlated at this location since only their proouct can be estimated. If the single thermo- couple were placed at any position, x, other than zero, then both k and pop can be simultaneously estimated provided the heating is continued long enough to cause the temperature to change at the sensor location. More complex mathematical models such as those describ- ing heatcxnnuct331 with temperature dependent thermal proper- ties or those describing heat conduction with change of phase do not reveal correlations between the parameters as simply as the example given above. This is because simple exact solutions are not known for these more general cases. Never- theless, numerical solutions can be utilized to find Optimum experiments in such cases. El? Sensitivity coefficients. Sensitivity coefficients can be utilized to aid in finding Optimum conditions. A sensitivity coefficient is defined as the derivative of the predicted value of a measured quantity with respect to the parameter pi, 3O arms») Sensitivity coefficient = 8g1(§,t) = BPi— [2.2.1] where f is a vector of controllable variables and P is a vector of parameters. T(P,£it) is the predicted response for the controllable variables set at f and time t. If for any fixed 6 and all values of t the relation Zijsgjgn) = 0 [2.2.2] 3 is satisfied and any one of the constant coefficients Aj is nonzero the sensitivity coefficients are said to be lin- early dependent. For cases where the sensitivity coeffici- ents are linearly dependent all the parameters cannot be found simultaneously<28). This is illustrated by differ- entiating equation [2.1.1] with respect to k to obtain F 1 S£(O,t) = Egg-Lil = - EQ-(t/‘Trpcpkrz [2.2.3] and with respect to Pop to obtain T _ aT§O,t2 _ _,39 i s Cp(O,t) .. 3P°p .. PCp(t/7rpcpk) [2.2.4] we see that for all values of t T '1‘ AlSk(O,t) + AZSPCP(O,t) = 0 [2.2.5] where A1 = k [2.2.5(a)] 31 and 12 = -pcp [2.2.5(b)] Any criteria for finding optimum experimental conditions should alert the experimenter to linear dependence of the sensitivity coefficients. 2.3, Criteria for design of Optimgm experiments. While there is much work in the literature on para- meter estimation, few researchers have considered the problem of experimental design. Seinfe1d(26) states that there may be cases where the parameters in a mathematical model can- not be estimated but does not consider the problem of de- signing experiments to avoid these cases. G.E.P. Box and co-workers seem to have been the first to realize the need for experimental design in parameter estimation. (29) Box and Lucas suggested maximizing the deter- minant of the matrix I‘ l‘. = 2.22. [2.3.1] as a criterioni The elements of the sensitivity matrix Q are _ T zij ’ Spi(2’6j't)£=’§ [2.5.1(a)_] is a least squares estimate of P, the vector of parameters. Pp) * Note the prime (') means transposition. . ... rug-wt." Dr. (I, [l’ A 32 If the errors in the measurements are independent and nor- mally distributed with zero mean and constant variance then the volume of the confidence region in parameter space is inversly prOportional to the square root of det 1:. By designing the experiment to maximize det I: the confidence region for P is minimized if these conditions are valid. These assumptions are not all valid in our case since the experimental errors are correlated rather than being independent. The nature of this correlation is not understood at this time. Other justification for maximiz- ing the determinant of.l: that is not dependent on these assumptions has been presented in (5). NO better criterion is available at this time. One justification of the criterion is to note that if any two or more sensitivity coefficients for a given experiment are linearly dependent as discussed in section 2.2, the criterion, det I1, is identically zero. This re- sults from the rows of the matrix ; being proportional(30). If det I: = O, the confidence region in parameter space is infinite whether or not all the assumptions are valid. Most models in heat transfer are nonlinear in the parameters. For such cases the sensitivity coefficients are dependent on the value of the parameters and thus an initial estimate of P is needed to design the best experiment for estimating P. It may seem strange that in order to use this Seheme one must have initial estimates Of the parameters that are to be found; however, any experimental design uses the experimenter's prior knowledge Of the subject. The better 33 his initial guess the better and more efficient the experi- mental design will be. The criterion can be made arbitrarily large by in- creasing the number Of observations and the range of the measured variable(s). Physically the experimenter is re- stricted to finite quantities, hence constraints must be imposed. For a large number of equally spaced measurements in time and a model that is linear in the dependent variable Beck<5> modified the criterion to include constraints of a fixed number of observations and a fixed temperature range. The resulting criterion is to maximize the determinant of the matrix 4; whose elements are n 1 - l m P T P [2.3.2] 15 is the maximum dimensionless duration of the transient experiment, n is the number of thermocouples, s; (xzyT) i and S? (x‘,77 are the sensitivity coefficients at location J :m‘ for parameters p1 and pj respectively. Zleax is the constraint of the maximum temperature rise occurringiin the experiment. Beck considered models that are linear in the de- pmnnient variables although the dependent variables may be runtlinear in the parameters. This can be a confusing but :hnportant distinction. For models that are linear in the chapendent variable the dimensionless sensitivity coefficients 34 are independent of temperature rise. Thus for models which are linear in the dependent variable it becomes unnecessary to consider temperature rise as a separate constraint. If the model is nonlinear in the dependent variable, the sen- sitivity coefficients are also nonlinear in the dependent variable. Consequently the choice of maximum temperature rise must be considered as a separate constraint. For non- linear problems there may exist an Optimum temperature rise. Heinekin, Tsuchiya, and Aris<3l) proposed a criterion for systems of ordinary differential equations that is similar to the criterion given by Beck. For the types of experiments conducted in heat trans- fer investigations, examples of controllable variables are thickness of sample, initial conditions, type of boundary condition (i.e., heat flux or temperature), maximum tempera- ture rise, and to a limited extent the functional form of the boundary conditions. These controllable variables must be adjusted to maximize the determinant of the matrix 4;. For a finite body the sample length, E, enters the optimization criterion through the dimensionless time, “2 If the thermal diffusivity, a, is not constant it can be evaluated at some reference temperature for convenience. Thus once the optimum value of Ti is known, the sample length can be chosen from other considerations and the maximum run time, t, can be found from equation [2.3.3] . For a semi-infinite body with only a single thermocouple 35 the thermocouple location, x, enters the criterion through the dimensionless time, 1' = a.t/x2 [2.5.4] After 15 is determined, the duration of the experiment and thermocouple location(s) can be chosen from other consider- ations. The procedure then is to fix all the controllable variables except ‘Th and calculate Z = det A [2.5.5] as a function of ‘Tm. The maximum temperature rise, [5T case of a step change in temperature at a boundary, this is max’ must be attained for each value of ‘Tm. For the automatically satisfied. For the case of a heat flux bound- ary condition; however, the heat flux must be adjusted for each value of ‘Tm so that the given stmax is attained in time ‘15. The maximum value of [S mam] = 230” [2.5.6] indicates the best value of 15 for that particular setting of the other controllable variables. This procedure is repeated for different settings of the other controllable variables until the best combination is found. The procedure is illustrated for several test cases below. Any Optimum experiment must be a compromise between theory and practice. Practical considerations play an 36 important role in the choice of the controllable variables. One of the most important of these is the choice of the functional form of the boundary conditions. Only those conditions which are readily attainable and measureable in the laboratory need be considered. Other practical con- siderations that should be taken into account when planning an experiment are the capabilities of the measuring equip- ment used. All measuring equipment introduces some error into the measurement whether in the form of random electron- ic noise or a biased error due to inaccurate calibration, for example. To minimize the effect Of measurement errors from noise on the calculated parameters, the overall temper- ature rise for the experiment should be kept as large as practical. This is introduced into the Optimization criter- ion later in the form of a "merit" function for the tempera- ture rise. When thermocouples are attached to a body there is always some question as to the exact location of the measuring junction. In Chapter III, the experimental tech- nique will be discussed and it can be seen that the possible uncertainty in thermocouple location is equal to the radius of the thermocouple wire used. This effect can be minimized by using relatively thick specimens and the smallest practi- cal thermocouple wire diameter. 2.4 Results for several sample cases. Semi-infinite case Several test cases were investigated to find the Optimum experimental conditions. The first case was that 37 of a semi-infinite body of liquid having constant thermal properties that freezes when subjected to a step decrease in temperature at the surface x = O. The initial fluid temperature is Ti and the surface temperature becomes To. The solution for the temperature distribution is(27) m (x t)-T T+ _ _§__’...__a.. .. .2..— .1. .Ts ' Ti * To ‘ em)” [(47);] [2'4'1] in the solid and T (x,t)-T (1-T+) t i '2? = ° = 1 - -—————;—T f [—9-] 3 Ti - To erfc[)\(a+) ] or c[ T :[ [2.4.2] in the liquid. The dimensionless position of the solid- liquid interface is given by 6"(7) = 53%)- : 2X7? [2.4.3] The constant ).is found as the root of ->\2 (afiikfii-ipe'xzdf erf()\) - T$erfc[7\(O-+)%] = wi'XL; [2.4.4] ‘where the five dimensionless properties are defined by T - 'r T; = ff—3efi [2.4.5] L. L; = P°p 8 Tf'To) [2.4.6] 1 = agt [2.4.7] :3 ll :1 m \ [:3 [2.4.8] k = tz/ks [2.4.9] The subscripts s and 1 refer to the solid and liquid re- apectively. The heat Of fusion is L and the fusion tem- perature is Tf- Figure 2.4.1 (a) shows the dimensionless temperatures Ti and TS, versus dimensionless time, T, for several values of the dimensionless heat of fusion, Lg, with the dimension- less fusion temperature, TE, held at 0.5. Note that the larger the heat of fusion, the longer it takes the tempera- ture at a given location to drOp below the fusion temperature. Figure 2.4.1 (b) shows the dimensionless temperatures ver- sus dimensiOnless time for several values of the dimension- less fusion temperature, TE, with the dimensionless heat of fusion equal to 1.0. Note that as T; approaches 1.0, there is a more pronounced effect of the phase change upon the temperature history. Sensitivity coefficients for all the properties are shown in Figure 2.4.2, for the common case of Figure 2.4.l(a) and (b). A number of observations can be made from an in- spection of this figure. For example, the curves for Si and s SidL have almost identical shape which means they are almost linearly dependent. Because this case of a semi-infinite body has these curves plotted in terms of a similarity variable 1’which includes both the independent variables x and t, an Observed linear dependence in Figure 2.4.2 39 means that the dependence is true for any nonzero x and all t's. If the body were finite, the sensitivities would not be a function of only one independent variable. Another observation is that the magnitude of'ST is small which c Pps means this experiment would be ineffective for estimating the heat capacity of the solid phase. Another Observation T T are not correlated indicating that this f is that SBLand S experiment would be suscessful in simultaneously determin- ing FL and Tf. If one calculates Ziopt for this experiment for all six prOperties as explained in section 2.3, the magnitude 22 is on the order of 10. indicating a very small chance of simultaneously estimating all six parameters. Three or four parameters could be estimated simultaneously, however. An examination of Figure 2.4.2 can tell which ones. As noted above ST P°ps upon the predicted temperatures. Hence one should not try to is always small and thus Pope has little effect estimate pops from such an experiment; if one does, then all the parameters calculated at the same time could be quite inaccurate. Also k8 and FL should not be simultaneously estimated due to near linear dependence. A group of three that could be estimated together is Tf, pL, and kl' By examining sensitivity curves such as Figure 2.4.2 one can learn a great deal in terms of possible Optimum experiments. One estimation problem is to find the density-heat of fusion product» FE” when the other properties are known. The controllable variables for this case are thermocouple location, maximum run time, initial temperature Ti,.boundary 4o .9310 'd c: w T; = 0.5 at m + g I- k = 1.0 .p m g 0.6 P 2' a) .. 8.0 [—4 U) 3 0'4" 2.0 H 8 + .a ' =o,5 3 LT 1.0 5,110.2 " Q 1 l l 1 j 0 2 4 6 8 10 Figure 2.4.1 (a) Dimensionless Time, 7' Dimensionless temperature of a freezing semi-infinite body subjected to a step decrease in temperature at the boundary for several value of the dimensionless heat of fusion, LT' 41 1.0 + m LT = 1.0 EH + v§ (1. = 1.0 Q 008 P k+ = 1.0 E4 a? E 0.5 - d H m p. a :3 + - m 0.4 h T? - 100 m 0075 m ":1 0.50 3. g 0.2 ’ 0.25 m S (3 ° 0 2 4 6 s 10 Dimensionless Time, 1’ Figure 2.4.l(b) Dimensionless temperature of a freezing semi-infinite body subjected to a step decrease in temperature at the boundary for several values of dimensionless fusion temperature, TT. 42 0.5]- E490 Q0 0.2 h- 0 P°pi ”Tf 0.1 h 6* Q 1 “—- U) 1 *- g 0.0 - .3 + “i + a, T = 0.5 g -o.1 T 0 of = 1.0 >. 15 k‘“ = 1.0 cu £3 3 L l l 4 '°°3o 1 2 5 4 Dimensionless Time,‘T Figure 2.4.2 Dimensionless sensitivity coefficients for a freezing liquid semi-infinite body with a step decrease in temperature at the surface. 43 temperature To’ and maximum temperature rise. In order to find an Optimum experiment for estimating pL, each of these controllable variables must be considered. Thermocouple location and maximum run time are both included in the dimensionless time given by equation [2.4.7] (for a single thermocouple). Hence maximizing Zwith respect to 7 con- siders both location and run time for a single thermocouple. The initial temperature and boundary temperature are included in the dimensionless parameters T; and L;. The range of T; is from zero to one. If T; is zero, i.e., the boundary temperature at x = O is equal to the fusion temperature, no freezing can occur and hence no information related to the heat of fusion could be contained in the data. Table 2.4.1 shows the value of ZSOpt with L$=l for several values of TE. For a single thermocouple, TE=1 (the initial temperature Ti equal to the fusion temperature Tf) seems to be a desirable experimental condition. The value of L; can be controlled by an appropriate choice of the quantity (Tf-To), or since T; will be set equal to one, (Ti—To). The best value of L; may be chosen by examining lsopt 4. versus LT’ we see that versus LE. Figure 2.4.3 shows Zlopt the Optimum experiment for a single thermocouple with T$=l is one for which L; is larger, or in other words, (Ti-TO) should be chosen as small as practical. At this point, one may decide just what the Optimum value of [STm x=(Ti-To) should be given some information a about the experimental equipment. In Chapter III the ex- perimental equipment will mediscussed. The maximum 44 —' + Table 2.4.1 Aopt for several values of TT with pl. the parameter to be estimated for the semi- infinite example. L; = 1 of’= 1 k+ = 1 TT 5 opt Topt 0.00 0.00000 ---- 0.25 0.00015 8.86 0.50 0.00512 5.84 0.75 0.01476 2.52 1.00 0.05709 1.91 45 0.08 0.06 Opt 0.04 0.02 0.0 0 ~ 2 +4 6 8 LT Figure 2.4.3 ZOpt versus dimensionless heat of fusion for a single thermocouple with (L the property to be estimated. 1.00? 0.75" M 0.50“ 0.25h 0.0 J l l l 0 100 200 300 400 AT max Figure 2.4.4 Merit function, M, versus ZkTmax 46 temperature rise for the equipment at the highest gain set- ting Of the DC amplifiers is about 370°F. Electronic noise introduced by various components is on the order of about 50F. Thus if our experiment were planned with a 10F tem- perature rise say, any information in the signal would be masked by the noise in the signal. To aid in selecting the optimum temperature rise, "(32) for the temperature we now define a "merit function rise based on the characteristics of the equipment used. This idea of a merit function for the experimental equip- ment is new in the field of optimum experiments for para- meter estimation. Since the noise introduced into the signal is not a function Of the temperature rise but re- mains constant, a linear function was selected. The value of the merit function at a szmax of 00F is chosen as zero and at a lleax of 370°F the value of the merit function is chosen as one. f the temperature rise goes beyond 370°F the output of the DC amplifier goes beyond the range of the analog to digital converter in the computer and a lower gain.must be selected; it is not convenient to do this in the middle of a test so the merit function is assigned a ‘value of zero for [BT greater than 3700F. Figure 2.4.4 max shows the merit function, M, versus ATmax' To select the Optimum temperature rise, the criteri- cni Of the maximum of the product of the merit function and 1aiopt is chosen and is shown in Figure 2.4.5. A [BTmax of aixn1t 350°F is the best using this criterion. If one wishes to calculate both the density-heat of 47 0.03r 0.02 MAW,c 0.01 0.0 Figure 2.4.5 Product of merit function, M, and zOpt versus ATm ax used to determine the Optimum ATmax for estimating QI" 48 fusion product and the fusion temperature using this ex- periment, Table 2.4.2 indicates that T$=l and large L; still provide the best experiment. The Optimum temperature rise was not computed for this case. Finite solid insulated at x=E The next example is that of a finite body of solid material that is insulated at the surface x=E. A constant heat flux, q, is applied at the surface x=O and the solid melts. Sensitivity coefficients for this case were calculated using the finite difference methods described in Appendices A and B. The dimensionless heat of fusion is given by L§=7584QL'EE7E, [2.4.10] EL and the dimensionless fusion temperature is given by T T T; = 3872: [2.4.11] Figures 2.4.6 through 2.4.8 show the dimensionless sen- sitivity coefficients for values of L;=l and T;=0.5 at locations x/E=0.0, 0.5, and 1.0 respectively. NO attempt ‘to find the Optimum experiment for this case was made. Visual examination of the sensitivity coefficients in Figures 2.4.6 through 2.4.8 shows. that all the curves Ihave distinct shapes at each location which means that they are linearly independent at least in pairs and thus pairs of parameters are not correlated. Closer examination reveals that the sensitivity coefficient SIT: has a significant 8 I Ir- 1‘ ,; ."Y." ‘1‘. luml I! e” SIP—J 49 Table 2.4.2 l>opt for several values of T; and L; with Tf and pL being the parameters to be estimated for the semi-infinite example. OP'==1.0 k+ = 1.0 TT LT zsopt prt 0.0 --- 0.0 --- 0.5 0.5 5.857x10‘5 1.60 0.5 1.0 1.310x10‘4 2.94 1.0 0.5 2.720x10"3 0.675 1.0 1.0 6.627x10'3 1.08 50 2.5" 2.0- ~14 a|a .920 .a 1°5T' t 9“, N 'F' 1.0- 9‘90 03 Dimensionless Sensitivity Coefficient, c -l.5 - < p” -2.0i— -2.5 1 l l 1 in I 0.0 1.0 2.0 3.0 Dimensionless Time , "C Figure 2.4.6 Sensitivity coefficients at x/E=0.0 for a melting-freezing, finite body with constant heat flux, q, at x=0 and insulated at x=E. 51 2.5— 2.0” H .4... are a 1.5“ t n'cr u 3.4 1.0- p13 k8 T f k 0.5— 1‘ 0.0 1 M J I Dimensionless Sensitivity Coefficient, -0.5- (L I 9° c ecps p8 (c 9 ps -1.0.. 1’" QL -1.5 '- (cps -2.0[— _2.5 J L l . J l _l 0.0 1.0 2.0 3.0 Dimensionless Time, Figure 2.4.7 Sensitivity Coefficients at x/E=0.5 for a melting-freezing, finite body with constant heat flux, q, at x=0 and insulated at x=E. 52 2.5r- 200 '- ...; sale. 200 p, 1.5e t p'cr u *4 l.0" ebfh T 4; 1‘4 a 0) H O ...; 44 4. 8 1 1 1 .4 :5 >3 :3 Ti b ...4 v I 7-3 c g e ps a: m U) Q) ’3 o 73 s -l.5 - 3 (Op; «4 Q L -200_ Q -205 1 l l L L J 0.0 . 1.0 2.0 3.0 Dimensionless Time, Figmme 2.4.8 Sensitivity coefficients at x/E=l.0 for a melting-freezing, finite body with constant heat flux, q, at x=0 and insulated at x=E. ..t , e.“ - an“, H. war—:4!- 53 magnitude for only a short time relative to the other para- meters and thus may be the least accurate of all the para- meters estimated from this experiment. Experience has shown that it is generally easier to estimate accurately a few parameters simultaneously than many. Thus we will concentrate on estimating the para- meters pL and Tf simultaneously. These two parameters are directly involved in model building and discrimination for the experiments analyzed herein. The other parameters will be measured using other eXperimental techniques which are described in Chapter III. T At the heated surface, x/E=0.0, both 30L and ST T f have nonzero magnitudes and distinct shapes with 33L finally approaching a constant value of about -l.0 and sgf going to zero after the entire body melts. This would be an excellent thermocouple location to estimate FL and Tf simultaneously for these experimental conditions. At x/E=O.5 both ST and pL 3% again have nonzero magnitudes and different shapes. f Both are somewhat smaller in magnitude for this location than at x/E=0.0. This location is good but not as good as the heated surface. At x/E=l.0 SEE and'Sg are again f not correlated but 33L has a very smallamagnitude until the entire body melts and then becomes very large for a short time. ng has a fairly large magnitude until the entire body melts and then becomes zero. Thus x=E is also a very good thermocouple location, if the experiment is allowed to run until some time after the entire body melts. Thus one can conclude that the heat of fusion FL, and the fusion temperature, Tf, can be estimated simultaneously 54 using these experimental conditions and a single thermocouple at many locations in the body. If more than one thermo- couple is used, the accuracy of the estimate can be improved. Figure 2.4.9 shows the dimensionless sensitivity co- efficient figm calculated in the same manner as the previous figures but with La=0.445 and T;=0.0. It is apparent from Figures 2.4.6 through 2.4.8 and Figure 2.4.9 that if one wanted to estimate only the heat of fusion, PL, using a constant heat flux at one boundary and the other boundary insulated, the most desirable location for a thermocouple would be at x=E, the insulated surface. The sensitivity co- efficient 33L has, by far, the largest magnitude at this location. An approximate solution to the freezing-melting problem is used to study this case in more detail. New approximate solution An approximate solution can sometimes give greater insight into a problem than a more accurate solution that is unwieldy. For this reason, an approximate solution is obtained using an integral energy equation method given by (33) Goodman The simple approximation of a straight line for the temperature distribution is used. In the approximate solution the body is initially a solid at the fusion tempera- ture, Tf. The solution is 11-53 - %-+ L” + [(L")2 + 21]* O§x$ em; . f = q q 7573 ‘13”; 0 €(T)éx_‘.E; 1‘1 1 [2.4.12] where 6(7? is the position of the liquid-solid interface [1 55 '1'41- — 1-21— x/E 8 IO ...-Op -03; ‘3 —O.6b ("6 ..J -0.‘ i- U. . 11 0'5 r-ma' -02... 04) 0.0 1 l 4 J 04) (>4 (>8 P2 P6 T nest/E2 Figure 2.4.9 Dimensionless sensitivity coefficient, 83L, versus dimensionless time,‘T, for a finite body with constant heat flux, q, at x/E=O and insulated at x/E=l. 56 and is given by + _ 6(7) __ __ + + 2 4‘; s (1) - a — Lq + [(Lq) + 21] [2.4.15] £5 ‘1' is the dimensionless time, 1' = -—2— . [2.4.14] L; is the dimensionless heat of fusion and is given by IDL + I‘q = Pop” (qE/E‘ ) [2°4'15] At the instant the entire body is melted, the tem- perature distribution is still linear; an exact solution for the problem of a body subjected to a constant heat flux on one side, insulated on the other side and with a linear tem- perature distribution at T = 1E is an T-T ‘ n 2 2 f 2 E [-1) -n 11’ (7-7“) nfix = 1’+ - - e n cos ; qE7k 2 'fl’n n2 ( E ) OéxéE, ‘T>‘TE [2.4.16] 13518 the time taken to melt the entire body and is obtained fixnn [2.4.12] by letting x=E and T-Tf=0; the result is t1) 1 = i— + La. [2.4.17] 'The approximate solution [2.4.12], is differentiated with respect to pL to find ST (x,t), 57 L. L” 05x5 em. .. + 3 q l + g; ‘7‘:1E L PL 7’ " qE7k 3]PL[ " o; €(T)£XéE, 7‘71: [2.4.18] For times greater than ‘TE the solution is found from [2.4.16] to be 2 2 T - S : - 11' (1-7 ) 1 fipL(x’t) ' ‘L; ' 2L; ('l)ne n E cos [nE x) The only thermocouple location considered here is the insulated surface. Only times greater than ‘TE are of interest because Sng’t) is zero until 7:71;. The maximum temperature rise in the liquid occurs at the heated surface =0 and can be shown to be AT O 2 2 327% = (T-T'E) + 2- 37223-14; e‘n 7' (7" E) [2.4.20] 11 n: Equations [2.4.19] and [2.4.20] are employed in calculating [l as described in section 2.3. Figure 2.4.10 shows A versus dimensionless time for several values of ATmax Remember that the maximum temperature rise,Z§T is an ex- max’ perimental constraint. For each value of 15 the heat flux, q, must be adjusted so that tsTmax is attained in time ‘Tm. To determine the Optimum temperature rise for this experiment the value of ZSopt i.e., the value of £5 at the optimum ‘Fm, is plotted versus lBTmax in Figure 2.4.11. 58 8— ATmax = 50 ] 6A- 4)- loop 21- 200[" 400 o r I I I I I 0 2 4 6 Dimensionless Time, 17111 Figure 2.4.10 ZS versus —rm for several values of liTmax for a melting, finite body of solid initially at the fusion temperature with (L the para- meter to be estimated. 59 Figure 2.4.11 Aopt versus ATmax for finite body of solid initially at the fusion temperature with {L the parameter to be estimated. 60 Figure 2.4.11 indicates that ATmag—o is the ideal condi- tion; however, arguments have already been presented to show that this is not the case in practice. To complete this case study, we refer to the merit function Of Figure 2.4.3. Again the optimum criterion is to be the maximum value of the M and ZSOpt product. The product M x [SOpt versus Zle ax is shown in Figure 2.4.12. From this figure the optimum temperature rise is found to be about 33°F. The Optimum dimensionless time,‘T is now fixed m,opt’ and can be determined by calculating ZSOpt for a [STmax of 33°F as in Figure 2.4.10. Given.7h’opt, the value of E can be selected and then q can be determined from equation [2.4.20]. It should be emphasized again that if it were not for the existance of a relation between the maximum temperature rise and the heat flux, such as equation [2.4.20], this optimization process would be more difficult. In the case Of a more accurate or more general solution such as a numerical solution, no such relation exists. The process thus becomes a more time consuming iterative type calcula- tion. For each value of Ti some other method of finding ‘the heat flux that will cause a maximim temperature rise AT costly preliminary procedure, the lessons learned from this max in time‘Tm must be found. Instead of attempting this simple model and visual examination of the sensitivity co- efficients will be used to design our experiment. Heat flux approximating_the actual case One other test case was examined. The actual form of 61 105 '- l.0 - .p a. o h- 14 I: 0.5 ' 0.0 - 1 1 l 1 O 100 A 200 o 300 400 . Tmax in F Figure 2.4.12 The product Mzopt versus ATmax used to determine the Optimum Mmax for estimating the parameter QL. 62 the heat flux expected from the experiment was used to calcu- late the sensitivity coefficients S£m(x,t) and 3; (x,t) for f a finite body with the surface x=E insulated. The heat flux used is that shown in Figure 2.4.13. The sensitivity coefficients igL(x,t) and 8% (x,t) versus time are shown for x/E = 0.0, 0.5, and 1.0 In Figures 2.4.14, 2.4.15, and 2.4.16 respectively. Examination of these figures again indicates that the heat of fusion, pL, and the fusion tem- perature, Tf, can be estimated simultaneously from this ex- periment since the sensitivity coefficients igL(x,t) and Si (x,t) are not correlated at any of the three locations shown. Table 2.4.3 gives the locations and times for the maxima and minima of the sensitivity coefficients 33L and T Tr infinite body with a step decrease in the surface tempera- S for all the cases that were computed. For the semi- ture the maxima occurs at the solid-liquid interface. The only zero values are at the surface x=0 which is what one would expect because the temperature is prescribed there and changes in the prOperties can have no effect on it. For the cases using finite bodies all the maxima occur'at either the insulated surface or at the heated sur- face and in most cases at times after the entire body has melted. Zero values of both 85L and 8,]; occur but not at the f sanue time. 33L is zero (actually not identically zero but relmitively small) for times less than Ty and 8% is zero for “ f tinues greater than 7%. figL being zero while 8% is not and f vice versa permits the simultaneous estimation Of pL and Tf. 63 60 o g H t Cal 1 t ea er- or me er +3 44 \\ a +3 an 53 3O - '3 Lid a. p m o a; 15_‘ OZ 1 1 1 0 75 150 Time in Seconds Figure 2.4.13 Experimental heat flux from OOpper heater-calorimeter and heat flux through copper lid versus time. 64 1000 — a 6.” ae[0 7.5+- [.34 v4\\ 94 It! My II 5.0r- H 849. a) S 2.5 - o .8 - T -H pi ' f 2: m M <8 0'0 1* 7r ‘\\='—? *1 >6 :1 5 9L k E a) ‘2-5 [— a Q) a) {D U) Q) :3 -s.o—- O H U) m (D .S . D 4.5%- -l0.0 1 1 ' ‘ .O 1.0 2.0 3.0 4.0 Dimensionless time, t: att/E‘2 Figure 2.4.14 Sensitivity coefficients SEL and 3; at f x=0.0 versus time for a finite body insulated at x=E and the heat flux of Figure 2.4.12 prescribed at x=0.0. 65 10.0 ' e167 7.5 - 204) .a e+\\ ... s u 5'0 I «4 etc. a) «4 pi‘ f O «"3 E 0.0 K171 l L l J c; 'L_ >, ]\c_____. .p , ml :‘3 e 3 -2.5 l- s Q) a) (0 3'3 '3 -500 [- O H O) c Q) 3 G -7-5 *- -10.0 #1 41 2,1 .sJ 0.0 1.0 2.0 3.0 4.0 Dimensionless time, 1' = «Lt/E2 Figure 2.4.15 Sensitivity coefficients SEL and s; at f x/E=O.5 versus time for a finite body insulated at x=E and the heat flux of Figure 2.4.12 prescribed at x=0.0. 10.0 P set“ ’0'» 7.5 +- 2.5 ‘ Dimensionless Sensitivity Coefficient, 66 -5.0 I- QL -7.5 ’ _l0.0 l l 1 J 0.0 1.0 2.0 3.0 4.0 Dimensionless time, t= ott/ 2 Figure 2.4.16 Sensitivity coefficients 82L and s; at f x/E=l.0 versus time for a finite body insulated at x=E and the heat flux of Figure 2.4.12 prescribed at x=0.0. 67 .aeon oaseso one eats oe cease ea case In degassed commons anonm noumnh on» gown: as made uneducaenosau on» as penance ma my. mooadoamnosac one as confines ma :9. Had zev ads MPV add SPA Ham o.H 0.H mun as covsdsnda .ouw ed ~H.¢.~ oedema no mean econ .seon ceases HAN Haw o.H e o o.ou+a .mee.ou+n mun es eoeoasesa ouw es wean econ assessoo .seoa ceased Add HHd Add 0.0 o.H o o m.on+a .Hu+a we coeoascsa so mean secs .seon ceases aux ouw assesses Ham 0.0 In. H60“ a». HOG“ assess Had aaaoaeasa ouw pd sedaqvmoo assesses onapauoasov soon oeaaamsauaaom M\N hr m\w hr M\w ur m\H m ohms «Wm owned a flow omhdd Hy bMUCH ammo [I musmaoaummoo haa>avamdom Ho «sands can saunas you H a nose» use wdoavdooq . am and Au 9m m.e.m manna 68 The sensitivity coefficient SFT1L has a relatively large magnitude at all locations in the body for times greater than ‘TM. The sensitivity coefficient sg has a f relatively large magnitude for all locations up until time T E. From the above considerations one can draw several conclusions as to the best experiment to estimate the para- meters pL and Tf simultaneously. 1. A finite body with one surface insulated and a .known heat flux at the other surface is desired. 2. The best locations for measurement are at the heated surface and the insulated surface. 3. The experiment should last somewhat longer than the time it takes to melt the entire body. 4. To estimate the parameter pL only,ATmax should be kept as small as possible within the limits of the experimental equipment. This may not be true if one wants to find both pL and Tf simultaneously. These conclusions were used to design an experiment to esti- mate the parameters in a melting—freezing model. The experi- ment will be described in Chapter III. CHAPTER III EXPERIMENTAL EQUIPMENT AND TECHNIQUE 3.1. The data acquisipion system. The data acquisition system of Michigan State University Thermal PrOperties Measurement Facility con- sists of three basic parts: Thermocouple sensors, a com- puter signal conditioner, and an IBM 1800 computer. The type and mounting technique of the thermocouple sensors is discussed in Section 3.2. The computer signal conditioner contains an electron- ic reference junction compensator and a DC amplifier for each thermocouple. The electronic reference junction com- pensators are made by Consolidated Ohmic Devices, Inc. model number EZT 213-A9. The function of the electronic reference junction compensator is to add to the voltage pro- duced by the measuring junction of each thermocouple a voltage that compensates for the actual reference junction temperature not being at the standard 32°F. Thus, the elec- tronic reference junctions are supposed to compensate for changes in ambient temperature and eliminate the need for bothersome ice baths or alternatively, calibrating before each use. It was found that they did not perform as ex— pected and calibration before each experiment was necessary. 69 70 The DC amplifiers are the Dana Laboratories In- corporated Model 3400. The purpose of the amplifiers, which have a maximum gain setting of 1000, is to boost the voltage produced by the thermocouples, which is in the millivolt range, to a level which the computer can read accurately. Some of the specifications for the Model 3400 are (54) frequency response 'DC to 100 Hz - I 0.01%, linearity AC to 2 kHz - 10.01%, noise at a gain of 1000 - 4 microvolts, input impedence - lO megohms. The low noise figure is highly desirable and the large imput impedence allows great free- dom in the choice of thermocouple wire size and lead length. There are nine separate reference junction-compensator-DC amplifier combinations allowing up to nine separate channels of data. The IBM 1800 computer is equipped with an analog to digital converter with multiplexer and an interval timer that allow the data from each thermocouple to be auto- matically discretized and recorded. The minimum time step between data points is about 50 milliseconds. Data is stored on a magnetic disk storage unit and the maximum number of data points is limited only by the available disk storage of the computer. A typical run may contain 2500 data points; without the computer, it would be a nearly impossible task to discretize and record this many data points. It was discovered that the reference junction com- pensators introduced 60 Hz noise to the signal to be ampli- fied. Two readings from each channel are taken 1/120 of a 71 second apart and then averaged to eliminate this 60 Hz noise. This was found to be a very effective method of eliminating the 60 Hz noise since the frequency of the signal from the thermocouple is much much less than 60 Hz. The data acquisi- tion system is calibrated on-line before each experiment. Four separate constant temperature sources are employed in the calibration, they are 8 Leeds and Northrup thermocouple checking and calibrating furnace Model 9009, boiling distilled water, a large copper block at room temperature and an ice bath. The furnace is set to a temperature of about 360°F and a Leeds and Northrup certified Platinum versus 10% Plati- num-Rhodium thermocouple with a Leeds and Northrup guarded six dial potentiometer Model 7556 are used to read the true temperature of the furnace. The true temperature of the boiling water is read using a Tagliabue certified mercury-in- glass thermometer with a range from 167°F to 221°F and grad- uations of 0.2OF. Similar Tagliabue certified mercury-in- glass thermometers with apprOpriate ranges were used to read the temperature of the copper block at room temperature and to check the temperature of the ice bath (32°F). Thermocouples made from the same spool as used in the specimens are also placed in the furnace, the boiling 'water, etc. These thermocouples are connected to a Leeds and Northrup rotary switch Model 8248 which allows one at a time to be connected to the computer signal conditioner. The computer reads each thermocouple ten times and stores the average voltage. When all four calibration points are read for all nine channels, the least squares quadratic, 72 T = A(i)v(i)2 + B(i)v(i) + 0(1) [5.1.1] is passed through the data. v(i) is the voltage read by the computer for channel 1 and T is the temperature read from one of the four standards. Thus the results of the calibration procedure is a set of coefficients A(i), B(i), and C(i) for each data channel. When the experiment is run, the voltages, v(i), are read and stored on the disk and equation [3.1.1] is used to convert them to temperatures at a later time. §.2 Egpipment and test specimens. The test specimens consist of short, solid cylinders of the "unknown" material three inches in diameter and of various heights. The cylinder is insulated on the peri- meter and at the bottom. Thermocouples are placed on the outer radial surface of the cylinder at various locations. Several thermocouples are located at the same height around the perimeter of the cylinder and connected to the com- puter signal conditioner in parallel to read an average temperature for that leVel. Thermocouples are attached to the surface of the cylinder in various ways depending on the type of material. For metallic specimens a small groove 0.010 inches wide and 0.010 inches deep is machined in the specimen at the desired point of attachment; number 30 gage wire (0.010 inches in diameter) is placed in the groove and the sides of the groove peened over to pinch the wire and hold it in place. A junction between the two dissimilar metals of the thermocouple wire is not formed 73 before the wire is attached to the specimen; the metal in the specimen becomes part of the actual measuring junction. This method eliminates the need for the usual adhesives which add mass and may possibly insulate the thermocouple from the specimen. For non-metallic specimens a junction is formed by electrically welding the thermocouple wire to form a "bead"; a shallow hole is drilled in the sample and the bead inserted and glued in the hole. Non-metalic sub- stances usually have much lower thermal conductivity than metallic substances so the adhesives usually will not cause problems by insulating the thermocouple since the conductivi- ty of the adhesive may be chosen to be near that Of the specimen. The specimen is mounted on the hydraulic cylinder in the loading frame shown in Figure 3.2.1. Mounted in the loading frame above the specimen of "unknown" material is a similar specimen of OFHC c0pper. The copper specimen or calorimeter is brought to some elevated temperature with the electric heater shown in Figure 3.2.1. When the desired temperature is reached, the bifurcated heater moves out from between the two specimens and the hydraulic cylinder brings the unknown and the standard into contact; the switch that controls the hydraulic cylinder also gives a command to the computer to initiate data aquisition. The sides of both the Specimen and the calorimeter are insulated and the heat flux from the copper is nearly uniform over the surface. Thus a one dimensional heat flow analysis can be utilized. For specimens which have a relatively high thermal 74 Copper Calorimeter Loading Frame Electric Heater Specimen Hydraulic Cylinder Figure 3.2.1 The experimental equipment. 75 conductivity the assumption of perfect insulation is reason- able. The heat losses from the sides of a high conductivity specimen such as used in this research are a small fraction of the heat flow through it. The applied heat flux over a surface is made as uniform as possible by carefully prepar- ing flat specimens and providing a uniform layer of silicon grease. A ”comb" can be made to provide uniform film thick- nesses of about 0.015 inch. One can choose either a temperature boundary con- dition or a heat flux boundary condition. If a temperature boundary condition is chosen, the data from the extreme thermocouples are used as the boundary temperatures. If; however, a heat flux boundary condition is required for parameter estimation, the procedure is more involved. The copper calorimeter whose thermal properties are well- known is instrumented with thermocouples for this case. We have found that the one dimensional heat conduction ‘equation describes the temperature distribution adequately. Thus the describing differential equation, the temperature history at several locations, the initial conditions, and the insulation boundary condition are all known for the copper calorimeter. The problem is now the inverse of the normal boundary Value problem; that is, the boundary con- dition at the heated surface is to be found using the given temperature history. Beck(35) has developed a successful technique for solving this problem and it is discussed in Appendix C. 76 3.3 The experiment. The experiment chosen to illustrate the model build- ing and discrimination techniques proposed in Chapter I was designed to take maximum advantage of the previously developed experimental techniques and equipment available. It was decided to illustrate the techniques for a melting material which is more complex behavior than linear heat conduction. A eutectic alloy Of 50% bismuth, 26.6% lead, 13.3% tin, and 10% cadmium was selected; this alloy melts or changes phase at 1600F which is well within the range of the equipment available. A nylon cup was designed to contain the liquid metal and hold the thermocduples in place. A % inch thick c0pper lid was pressed into the top of the nylon cup to keep the liquid metal from overflowing when the OOpper heater was pressed into place, see Figure 3.3.1. A nylon plug with pipe threads was used to seal the filler hole in the side Of the nylon cup. The thermocouple wire was 26 gage iron versus constantan. These thermocouples gave a voltage of about 10 millivolts for a temperature of about 370°F; for a gain of 1000 the computer signal conditioner will then supply a maximum of 10 volts to the computer. The IBM 1800 computer at MSU can digitize signals only between t 10 volts. The wire size chosen was 26 gage or 0.0159 inches in diameter. Two 0.0160 inch diameter holes were drilled through the side of the nylon cup at the location of each measuring junction; the iron wire was inserted in one hole and the constantan wire 77 .moamdoooahonp MO noawsooa was woman ma cad Monaco weakens condemns use dedhd no zma> soapoom Heapasm H.m.n madman muoaedooq odmsoooauona we no noahz Hm dash / 7r-|t O O 72.6 e 8.90 . o o \‘\ Ibmmwé Huh 0 o o V \ redo. % O O / . f C 1 N \ / menace hoadd ohfivdhonaoa son Lr 78 inserted in the other and both were epoxied in place as shown in Figure 3.3.1. The measuring junctions were com- pleted when the nylon cup was filled with liquid metal. The nylon cup was insulated on the sides with transite and fastened to an insulating base, also made Of transite, which allowed it to be mounted on the hydraulic cylinder. The conclusions obtained in Chapter II were utilized in designing the experiment to estimate the parameters in the melting-freezing metal alloy. A finite body with a known heat flux at one surface and insulated at the other surface is desired. The insulation boundary condition is approximated as closely as possible by the low conductivity material used for the cup. (Nylon has a thermal conductivity of about 0.14 Btu/hr/ftoF). The heat flux from the copper calorimeter-heater can be calculated as explained in Ap- pendix 0. As shown on Figure 3.3.1, there are thermocouples as close as physically possible to the insulated and heated surfaces which are the best locations for measurement. The experiment should run for some time longer than it takes to melt the entire specimen. The actual run time was determin- ed from this criterion and the thickness of the sample. The sample thickness was chosen large enough to keep the effect of the inaccurately known thermocouple locations to a mini- mum but small enoughsm that enough energy could be stored in the copper heater to melt the entire specimen. The thick- ness was chosen as 0.75 inches and the experimental run time was then calculated from the optimum conditions described in Chapter II to be about 150 seconds. 79 The computer program for calculating temperature distributions in the case of melting-freezing problems with a heat flux boundary condition requires the front boundary condition to be applied to the melting or freezing material, thus the heat flux leaving the copper lid must be calculated. The copper lid cannot be considered as part of the c0pper heater-calorimeter in the heat flux calculation because contact resistance between the two causes a discontinuity in the temperature distribution. To avoid this problem, the heat flux from the copper heater- calorimeter was computed separately using the method in Appendix C. The method in Appendix C was also usedtto compute the heat flux leaving the c0pper lid. The tempera- ture history in the copper lid was obtained from a thermo- couple at the surface nearest the cOpper heater-calorimeter and the boundary condition at that surface is assumed to be the heat flux that leaves the copper heater-calorimeter. Experience has shown that the larger the number of properties that are estimated simultaneously, the greater is the probability that there is correlation between them. This results in inaccuracy. To reduce this possibility, a second specimen was constructed to be used in estimating the thermal conductivity and density-specific heat product of the solid phase. A three inch diameter by one inch high cylinder was cast of the low-temperature alloy. It was machined smooth and thermocouples attached as described in Section 3.2. [Beck(5) recommends that thermocouples be located at both the 80 heated and insulated surfaces for determining constant thermal conductivity and heat capacity simultaneously. Three thermocouple locations were used with this specimen, they were 0.0625 and 0.125 inches from the heated surface and 0.0625 inches from the insulated surface. Prior to using this sample it was soaked in a refrigerator so the initial temperature was about 45°F. This allowed a greater temperature rise during the experiment before approaching the melting temperature which was avoided in this test. A separate experiment was also run using the nylon cup specimen to determine the thermal conductivity and density-specific heat product of the liquid phase. The sample was heated to about 180°F initially and the c0pper calorimeter heated with the electric heater to 400°F. Since heating Of the sample is from above, natural convect- ion effects in the liquid are negligible. One test was run with the c0pper calorimeter at a lower temperature than the liquid and the results showed natural convection causes a 30%>higher "effective" thermal conductivity to be calculated. The final run was made to determine the fusion tem- [perature and the density-heat of fusion product (or for the melting-over-a-range-ofhtemperature model the range and the height Of the "Spike" in the density-specific heat product) and to obtain data against which to compare all the pro— [posed models. The nylon cup specimen was heated from above fronlan.initial temperature of 75°F to a maximum temperature of almost 3OOOF. CHAPTER IV TEST OF THE PROPOSED METHOD 4.1 Introduction In this chapter the methods for model building and discrimination proposed in Chapter I, the results of the sample cases on Optimum experiments for parameter estima- tion in Chapter II, and the experimental equipment des- cribed in Chapter III are used to illustrate the model building and discrimination procedure. Even though we know the sample material is a low melting point alloy, we initially assume ignorace of this fact. Data used to test the model building-discrimina- tion procedure was obtained as described in Chapter III. Both the heat flux from the OOpper heater and the heat flux into the specimen (through the copper lid) are shown in Figure 2.4.13. These were calculated by the method out- lined in Appendix C. 5.2 The model building process As a first step in the model building process, a simple constant thermal prOperty model was fitted to the data. This corresponds to Block 1. of Figure 1.7.1. Beck(5) has shown that the optimum experiment for a con- stant thermal property heat conduction model has one surface 81 82 heated and the other insulated with thermocouples at the heated and insulated surfaces. These conditions are met by the experiment described in Chapter III. This satisfies Block 2. of Figure 1.7.1. The experiment was performed and the parameters estimated as in Blocks 3. and 4. of Figure 1.7.1 respectively. The property values that minimized the sum of squares function (see Appendix B) are 5.82 Btu/ hr/ft/OF for thermal conductivity and 90.37 Btu/ftB/OF for heat capacity. We now proceed to Block 5. the examination of the model and the residuals. The data and the model shown versus time in Figure 4.2.1. From this plot one can see the lack of fit for all thermocouples away from the heated surface. This plot clearly shows that the model is inadequate but it does not present the information on the deficiencies in the model in a form that is easy to interpret. Figure 4.2.2 shows the residuals, the difference be- tween the experimental and calculated temperatures, plotted against time. Note that except for the first thermocouple (x=O) positive peaks in the residuals for each thermocouple occur in chronological order moving from the heated surface toward the insulated surface. This may indicate a front moving through the body causing large errors in the pre- dicted response. Figure 4.2.3 shows the residuals for each thermo- couple plotted against the measured temperature of that thermocouple. This plot shows the residuals below a tem- perature of about 160°F to be nearly identical in both «I‘ll-ll IL “lid-1|- ] F.’ . .- 83 300-0 - thermocouph no. P 225-0» . t‘- " 9 . ‘ 3 ' . 0 "@J m '- . ' , c5 . . . ' ............ t—ISO'Of- , .' , -':Z.. -- (l60°F pcp(T) = 12.24; Tl60°F ] Btu/ftB/OF Proceed to Block 5. of Figure 1.7.1; Figure 4.2.4 shows the model and the data versus time. The prediction °F TEMPERATURE IN 3000 '- 225-0 *- ISO-0 *- 75-0 " / 87 model: step change in thermal properties thermoc0uple locations km {12-24 r< l60 0?} a,” . ~ o 000" ' —— -- I = 0 IVS" . O 627 urea F hrtt’F 3 ”-0350. c T). 2"‘8 1‘ '60 .f a.“ o __ ‘ S o 525” 39‘ I820 r- I60 '7 n3 or 0 - 1 = o 700' E 'I 0‘ 725' 00 TIME IN SECONDS Figure 4.2.4 Step change in thermal prOperties model and data versus time. if." *7; 1' I." ~ q 88 from this model is closer to the data than was the prediction for the constant thermal property model. However, this model still lacks the characteristic shape of the data especially near the insulated surface. This lack of characteristic shape indicates that this model is incapable of accurately describing this heat transfer phenomenon. Figure 4.2.5 shows the residuals for each thermo- couple versus time. This figure is harder to interpret than the previous examples. With the exception of thermo- couple number one at the heated surface, this model is an improvement over the constant thermal property model except at later times near the insulated surface. The peaks in the residuals for each thermocouple still seem to occur in chronological order but are negative instead of positive. Figure 4.2.6 shows the residuals for each thermo- couple versus the measured temperature for that thermo- couple. Note that the peak error for each thermocouple still occurs at approximately the same temperature, near l6OOF. Moving on to Block 6. of Figure 1.7.1; the modification made to the constant thermal property model is insufficient to describe the heat transfer process. There is still some unmodeled phenomenon occurring. The calculated heat capa- city for temperatures greater than l60°F is 182.0 Btu/ftB/OF. This large value also leads to rejecting this model because it is several times higher than that for any known substance. With the information gained from the models proposed this far, we move on to block 7. of Figure 1.7.1 and pro- pose two additional models. Both Figures 4.2.1 and 4.2.4, the data and the predicted responses versus time for both 89 50-0F _ model: step change In thermal properties '— ATHMHZ-l'tt 25-0 9 0 3 h- t-0 o ./ I9 I’ 0 h” I / 0‘0m1 IYITf 1T1 rrflr .J ‘\ g 9 0 9 o :7: ' \‘ a F V \ P a -25-o- © 0 .. thermocouple no. 6 ~ 0 P “50.0 I l 1 l L l l I l l J l l l I L l l .1 I 00 37-5 75-0 "25 L300 TIME IN SECONDS Figure 4.2.5 Residuals for the step change in thermal prOperty model versus time. 90 50-0- P __ model: step change in thermal propertles 250- 0 9 3 .3 I 9 I8 I- / .- O'C - I T T I T l T RESIDUAL, T T? I 0 -250. thermocouple no. 0 l- 9 _50.0 J J l l L A L 4 L l l A L l l L L J l J 00 75-0 ISO-0 225-0 300-0 MEASURED TEMPERATURE IN 'F Figure 4.2.6 Residuals for the step change in thermal properties model versus measured temperature. 91 the above models, seem to indicate that energy is being absorbed and not accounted for by either model. Figure 4.2.6 indicates that the unaccounted for energy is being absorbed in a small temperature range around l6OOF. In order to account for this absorbed energy, two different models were proposed; an isothermal change of phase model and a melting over a range of temperature model. The melting over a range of temperature model was simulated by a triangular pulse in the heat capacity. The solution to temperature variable thermal property problems is discussed in Appendix A. Details of the pulse appear on Figure 4.2.7. The parameters kl, the thermal conductivi- ty of the solid, and pcpl, the heat capacity of the solid, were evaluated from the solid, non-melting specimen des- cribed in Chapter III. The parameters k2, the thermal conductivity of the liquid, and pop}, the heat capacity of the liquid, were evaluated from a separate experiment using the nylon cup specimen with the sample initially all liquid as described in Chapter III. The experiment perform- ed in Chapter III appears to satisfactorily cover the range of controllable variables so that Block 2. of Figure 1.7.1 can be skipped. The height of the triangular pulse fbcpz, the range of melting and the location of the apex of the pulse, Tf, were calculated from the data shown in Figure 4.2.8. Values for all the parameters obtained are shown on Figure 4.2.8. The numerical solution of isothermal melting- freezing problems is also discussed in Appendix A. The e°p2 -— —--——————-- I"> Energy // Absorbed +3 H 3 ° I I; I o I "3 I I I Qcpl : 1 T1 ’2. T3 Temperature Figure 4.2.? Energy absorbed by melting over a range of temperature model. 93 SOO'OF ....... l' . . . r- thermocouple no. P 2250*- r LL I- 0 z I 3:1 l- D 2 I50'OI- a: g P -- model :4 r- - model: meltlnq over a range of temperature ---- data I. . j range I I. — | thermoc0uple locations 750 o-;-oooo ‘2 --—f— e—x-o-ns ‘1' "F 0 —- n . 0350 P om Btu r 0 - x . o- 525 III - I2-24m I12 - 6-27W e — g c 0.700 .. u E = 0.725 I -4 I - . Q°pl 24 e793; ecu: 35325. .F Q‘ns'" «$9.; *' range - 6-0 W t, - Ieo-o °r 0c 1 J J LJILL J L J l I L l l J J l J 1 l J 00 37-5 75-0 "25 ISO-O TIME IN SECONDS Figure 4.2.8 Melting over a range of temperature model and the data versus time. 94 parameters ks_and facps were evaluated from the all solid, non-melting specimen. The parameters kl. and Pep; were evaluated from a separate experiment using the nylon cup specimen initially all liquid as described previously in Chapter III. Optimum experiments for this case are dis- cussed in Chapter II; this satisfies Block 2. of Figure 1.7.1. The parameters Tf and pL were calculated using the data shown in Figure 4.2.8. If one calculates the energy absorbed by finding the area under the triangular pulse as shown in Figure 4.2.7., it should be approximately equal to the heat of fusion in the isothermal melting-freezing model. The area under the Btu lbm heat of fusion was estimated from the isothermal melting- freezing model to be 10,625 §%fi; thus they agree to within 1%. triangular pulse was calculated to be 10,683 and the Figure 4.2.8 shows the temperature plotted against time. Note that there are oscillations in the calculated temperature; this is due to the sensitivity of the numerical scheme to large changes in the physical properties. The oscillations damp out with time indicating that numerical method remains stable. Proceed now to Block 5. for the melting over a range of temperature model. The improvment in fit for the thermocouples near the insulated surface is apparent from an inspection of Figure 4.2.9 which shows the residuals versus time for all five thermocouples. Thermocouple number one which is located at the 50-0 - 25-0 " Tap" Tcal RESIDUAL, —2eo» - 50-0 00 95 thermocouple no. G model: melting over a range of temperature Mn“,- 6-785 11J1111J1L_11111111l 37-5 7.5-0 Ila-5 ISO-O TIhE IN SECODDS Figure 4.2.9 Residuals for the melting over a range of temperature model versus time. 96 heated surface has the largest error at the early times. This same residual pattern in thermocouple number 1 was also observed in all previous figures discussed. This pattern has also been noted in many previous experiments conducted at the Michigan State University Thermal Properties Measurement Facility with constant thermal property models. The most probable explanation is non- uniform heat flux across the surface causing three-dimen- sional effects to influence the measured temperature near this surface. Silicone grease applied carefully and uni- formly can help to reduce these effects. Figure 4.2.10 shows the residuals versus measured temperature. We note that the residuals for thermocouples 3,15 and 5 show similar characteristics when the material is completely solid, i.e., their sign is negative with a peak in amplitude just before melting begins. After melting these thermocouples exhibit residuals that have the same sign and increasing amplitude. In the melting range we note small spikes in the residuals indicating that the model does not completely describe the complex behavior in the melting range. For the isothermal melting-freezing model, Figure 4.2.11 shows the data and the model versus time. This model is also a great improvment over the first two. We note that there is no oscillation in this model as in the melting over a range of temperature model. Figures 4.2.12 and 4.2.13 show the residuals versus time and measured tempera- ture for the isothermal melting-freezing model respectively. 97 50-0r- l- l- model: meltlng over a range at temperature I- 2.5-Ol- r- r- E F ’— ' b 9 I'" 00 T t T r J l' 4 D D 9 m m b a: y- -25-0l- r- » thermocouple no. 0 p _50. 4 L 1 1 1 L 1 J 1 1 4 1 1 1 1 1 1 l 1 1 0-0 75-0 ISO-0 225-0 300-0 ' MEASURED TEMPERATURE IN °F Figure 4.2.10 Residuals for the melting over a range of temperature model versus measured temperature. 98 300-0- 225-0I- .L I 0 g I' U P r: 3 P ’— <1 {5 ISOOI- g . .— model: meltlng- freezing (isothermal) " thermOCouplu locations , . Btu . Btu I s - 0000 750/ "s '2“ hrtt’F "I'Gflmrvr 2 «wows . ‘ ' . a.“ ‘ . a! 5 - I '- O 550 P QC: 2‘ 48—f'3°F QC‘ 32 62 "3.F o _ l , o 525 Btu 0 — u - 0.700 T = I60-0 'F I. . Ioees-o— . - I- , Q "3 E o 725 r- 00 L I l I. L I. I I l I l I l J l I J J I J 0-0 37-5 75-0 "25 ISO-0 TIME IN SEC DVDS Figure 4.2.11 Isothermal melting-freezing model and the data versus time. 99 50'0F i- in 25-0 - thermocoupie no. 9 3 b I! .. 0-0 . .J < 3 9 U') u: 0‘. G) 0 -250" model: melting - freezing (Isothermal) . 0 ' ATflm- 7-03l _50.c 1 j 1 J LL L I L 1 J. I I I I I L I J I J 00 37-5 75-0 IIZ-S ISO-0 TIRE IN SECONDS Figure 4.2.12 Residuals for the isothermal melting- freezing model versus time. 100 50-0 [ model: melting-freezing (isothermal) 25-0 - - cal rep 0 0 d ..I' <1 3 O )- (7) Lu _ a: - 25-0- thermocouple no. b r- p _ 500 I 1 l I JLI I J l L L I I I I L J L I I 0 0 750 I500 225-0 300-0 MEASURED TEMPERATURE IN “E" Figure 4.2.13 Residuals for the isothermal melting- freezing model versus measured temperature. 101 The residuals are very similar in character to those for the melting over a range of temperature model. Both the melting over a range of temperature and isothermal melting—freezing models appear to be inadequate after the entire body has melted; the residuals for all thermocouples versus time have a negative SIOpe. This in- dicates that the calculated temperatures rise more rapidly than do the measured temperatures. We can conclude that the heat capacity for the liquid phase used in the models is too small. In fact, our early attempts to estimate the thermal conductivity and heat capacity in the liquid phase and the heat of fusion and fusion temperature simultaneously gave values for the heat capacity of the liquid which seemed erroneously large; for this reason it was decided to estimate the properties of the liquid phase from a separate experiment. It was possible to run only one experiment to evaluate the properties of the liquid phase because the specimen was accidentally overheated and the nylon cup destroyed. The cause of the discrepency between the liquid phase properties evaluated from the spearate all liquid experiment and those obtained from the experimental data shown of Figure 4.2.8 is not understood. It was decided to accept the prOperties calculated from the all liquid experiment because they were closer to published data(36). We proceed to Block 6. for both models. Examination of Figures 4.2.8 and 4.2.11 shows that both models have the same characteristic shape as the data they are supposed to 102 represent. That is, both models follow the data reasonably well and are able to absorb energy at nearly a constant temperature as evidenced by the thermocouples near the insulated surface. Since both these models are plausible, we now proceed to the discrimination stage, Block 8. in Figure 1.7.1. 4.3 Discrimination between rival mathematical models To discriminate by the likelihood method described in Chapter I, an estimate of experimental error is needed. One statistically proper way to obtain an estimate of vari- ance is to perform replicate runs. A true replicate run is not just a repeat of the experiment but must be made with completely different specimens and the equipment must also be recalibrated. The experiment is repeated and the differ- ences in the measurements of the two or more replicate runs are then used to estimate the variance. This can be an expensive task but it is usually recommended. Another way to estimate eXperimental error in temperature measurements is to use an estimate based on previous experience with the experimental equipment. This latter method is the procedure used in this work. During the course of checking out the experimental equipment,several tests were run to determine the variation in temperature indicated by each of the nine thermocouple- reference junction-data channel combinations. The process <20nsisted of simply having the computer read the data signal lO} continuously for an arbitrary length of time while all nine thermocouples were imbedded in an isothermal block of alumin- um. Data from these tests showed that the variance between readings on the same data channel was never more than 0.50F and between channels was usually not more than 50F. With this data as a basis, 50F was chosen as an estimate of I experimental error. It should be noted that the test used as a basis for I this estimate of experimental error includes only the error L introduced by the equipment used to amplify and read the voltage produced by the thermocouples. Errors caused by thermocouple wire variation, inaccurately known thermocouple location, heat losses from the sides and bottom of the specimen, expansion or contraction of the material during melting, and other similar errors are not accounted for in this test. However, these conditions tend to cause biased errors; that is, they tend to be one-sided for a given test. One purpose of replicate runs is to eliminate these biased errors. With a single experiment one has no way to discern these biased errors other than using several thermocouples- at each position relative to the heated surface. For the purpose of comparing likelihood functions we will follow the example set by Reilly(22) and multiply the likelihood functions for each model by an appropriate, common constant to make their magnitudes more tractable. Table 4.5.1 shows the likelihoods for each model that have been normalized by dividing by the likelihood of the step change in thermal prOperty model. The table of 104 Table 4.3.1 Comparison of the likelihood functions of the four models. Model mrms L(model/data) W ggggzfifl; thermal 18.70 0.0009 0.0235 giggmgiagggpigties 12.74 0.0389 1.000 ggltigge::§firz range 6.785 0.3982 10.23 Isothermal melting- 7.031 0.3721 9.559 freezing 105 likelihoods shows that the constant thermal property model is an implausible choice as the magnitude of its normalized likelihood is very small. The normalized likelihoods for both the melting over a range of temperature model and the isothermal melting-freezing model show a marked preference for either of these models over the step change in thermal properties model. The normalized likelihoods for the melt- ing over a range of temperature model and the isothermal melting-freezing models show that no discrimination is possi- ble between them, in other words the two models describe the data equally well. This is not an unreasonable conclusion since the isothermal melting-freezing model is really a special case of the melting over a range of temperature model, i.e., the range is zero for the isothermal model. Note that the melting range for the melting over a range of temperature model is only 60F; the estimate of experiment- al error was 50F, thus it is possible that the true melting range could be near zero. For this case both of the two melting-freezing models are capable of describing this com- plex heat transfer process. Practical considerations lead one to choose the isothermal melting-freezing model because it contains one less parameter and the numerical solution does not oscillate as does the melting over a range of temperature solution. CHAPTER V CONCLUSIONS This thesis considers one of the fundamental problems in science - how to build improved mathematical models from observations. This has been given the name of model build- ing and uses concepts of discrimination. A procedure to help an experimenter "build" an adequate mathematical model for cases where there are large amounts of transient data available at relatively low cost has been proposed. The proposed procedure for model building has been illustrated by finding several rival mathematical models for a complex heat transfer process. The specific heat transfer process chosen to illustrate the procedure was a finite, one-dimensional, solid, melting body. A heat flux boundary condition was applied to one surface of the body while the other surface was insulated. Temperature data was obtained by actual experiment. Through the example it was demonstrated that the deficiencies in a mathematical model can be brought to light by examining the model and the residuals. Knowledge of the deficiencies in a model is the information the experimenter needs to pr0pose an improved model. A method to help an experimenter discriminate between several rival mathematical models has also been incorporated in the model building procedure; it is illustrated using the 106 107 same experimental data. The discrimination problem is a fundamental and difficult one. Some analytical work has been done on this problem as indicated by the literature review. There is a need, however, of more discrimination studies that attempt to extend and apply the concepts to actual experimental data. This thesis is intended to help satisfy this lack particularly related to heat transfer. It is shown that not all discrimination procedures can be applied to heat transfer problems due to the expense of running each test and the dynamic nature of the tests. The method proposed for model building including discrimination was successful for the specific case of the melting low temperature alloy. The best mathematical model for this case is the isothermal melting-freezing model. The method presented could be used in other situations to find a mathematical model for a process; that is, under conditions where: 1. Large amounts of data can be obtained easily; e.g., where there is a transient experiment with thermocouples at several positions in the body. 2. The optimum conditions for parameter estimation can be found; i.e., the paramaters can be estimated. 3. A solution to the mathematical model can be found; this may mean that a digital computer is necessary. 108 The methods presented must be tempered with practical considerations and common sense. If this rule is observed the methods presented here can be useful tools to the ex- perimenter who is seeking a mathematical model for a process. BIBLIOGRAPHY 10. 11. LIST OF REFERENCES Box, G.E.P., Hunter, W.G., "An Experimental Study of Physical Mechanisms," Technometrics, VII (1965), 23-42 Pfahl, R.C., Mitchel, B.J., "Simultaneous Measurement of Six Thermal Properties of a Charring Plastic." Int. J. Heat and Mass Transfer, XIII (1969). 275-286 Box, G.E.P., "Fitting of Empirical Data," Mathematics Research Center, United States Army, The University of Wisconsin, MRC Technical Summary Report #151 (May 1960) Box, G.E.P., Coutie, G.A., "Application of Digital Computers in the Exploration of Functional Relation- ships," Proceedings of the Institugggn of Electrical Engineers, CIII (1956), 100 Beck, J.V., "The Optimum Analytical Design of Transient Experiments for Simultaneous Determinations of Thermal Conductivity and Specific Heat," PhD Thesis, Michigan State University (1964) Beck, J.V., Dhanak, A.M., "Simultaneous Determinations of Thermal Conductivity and Specific Heat," ASME paper No. 65-HT-l4 (1965) Pfahl, R.C., "An Experimental and Analytical Investi— gation of the Internal Mechanisms of Ablative Heat Transfer in Charring Cork," PhD Thesis, Cornell University (1965) Pfahl, R.C., "Nonlinear Least Squares: A Method for Simultaneous Thermal Property Determination in Ablating Polymgric Materials," J. of Applied Polymer Science, X 1966 Box, G.E.P., Hunter, W.G., "A Useful Method for Model Building," Technometrics, IV (1962), 301-318 Cochran, W.G., Cox, G.M., Experimental Designs, John Wiley & Sons, Inc. (1950), 122 Anscombe, F.F., "Examination of Residuals," Proceedings of the Fourth Berkelenyympogium on Mathematical Statistics and Probability, 1(1961), 1-37 109 12. 13. 14. 15. 16. l7. 18. 190 20. 21. 22. 23. 24. 25. 26. 110 Anscombe, F.J., Tukey, J.W., "The Examination and Analysis of Residuals," Technometrics, V (1963), 141-160 Draper, h., Smith, H., Applied Regression Analysis, John Wiley & Sons, Inc. (1966) Hunter, W.G., Reiner, A.M., "Designs for Discrimina- tion Between Two Rival Models," Technometrics, VII (1965), 307-323 Roth, P.M., "Design of Experiments for Discrimination ?mong Rival Models," PhD Thesis, Princeton University 1966) Box, G.E.P., Hill, W.J., "Discrimination Among Mech- anistic Models," Technometrics, IX (1967). 57-71 Hill, W.J., Hunter, W.G., "A Note on Designs for Model Discrimination: Variance Unknown Case," Technometrics, XI (1969). 396 Meeter, D., Pirie, W., Blot, W., "A Comparison of Two Model Discrimination Criteria," Technometrics, XII (1970). 457-470 Chernoff, H., "Sequential Designs of Experiments," Ann. Math. Statistics, XXX (1959), 755-770 Hill, W.J., Hunter, W.G., Wichern, "A Joint Design Criterion for the Dual Problem of Model.Discrimination and Parameter Estimation," Technometrics, X (1968), 145 Atkinson, A.C., "A Test for Discrimination Between Models," Biometrika, MVI (1969), 337 Reilly, P.M., "Statistical Methods in Model Discrimina- tion," Third Canadian Symposium on Catalysis, (October (1969) Hoel, P.G.,"Introduction to Mathematical Statistics," John Wiley & Sons, Inc. (19627: 57 Beck, J.V., Parameter Estimation in Engineering and and Science, To be published. Kmenta, J., Elements of Egpnometrics, The MacMillan Co., New York, N. Y., (1971). Seinfeld, J.H., "Identification of Parameters in Partial Differential Equations," Chemical Engineering Science, XXIV (1969). 65-74 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 111 Carslaw, H.S., Jeager, J.C., Conduction of Heat in Solids, Oxford University Press (1959), 75 Beck, J.V., "Analytical Determination of High Tem- perature Thermal Properties Using Plasms.Arcs," Pre- sented to the Eighth Conference on Thermal Conductivity, Purdue University, October 1968 Box, G.E.P., Lucas, H.L., "Design of Ex eriments in Nonlinear Situations," Biometrics, MVI {1959), 77 Stoll, R.R., Linear Algebra and Matrix Theory, McGraw- Hill Book Co., (1959), 89 Heineken, Tsuchiya, Aris, "On the Accuracy of Deter- mining Rate Constants in Enzimatic Reactions," Mathematical Biosciences, I (1967). 115-141 Strodtman, C.L., "A Design Optimization Technique Applied to a Squeeze Film, Gas, Journal Bearing," PhD Thesis, Michigan State University, (1968) Goodman, T.R., Advances in Heat Transfer, Academic Press, I (1964), 51 Anon., Instruction Manual for Dana Series 3400 Data Amplifiers, Publication Number DL 370, (1964) Beck, J.V., "Nonlinear Estimation Applied to the Non- linear Inverse Heat Conduction Problem", Int. J. of Heat and Mass Transfer, XIII (1970). 703-716 Anon., "Cerro Alloy Physical Data and Applications", Cerro C0pper and Brass Co., Bellfonte, Pennsylvania Smith, G.D., Numerical Solution of Partial Differ- ential Equations , Oxford University Press, (1965), 20 Murray, W. D., Landis, F., "Numerical and Machine Solutions of Transient Heat-Conduction Problems in- Yolving Melting or Freezing", ASME Paper No. 58-HT-16, 1958) Heitz,W.L., Westwater, J.W., "An Extension of the Numerical Method for Melting and Freezin Problems", Int. J. of Heat and Mass Transfer, XIII ?l970), 1371 Deutsch, R., Estimation Theopy , Prentice-Hall Inc., (1965). 60 APPEN DI CES APPENDIX A FINITE DIFFERENCE METHOD FOR CALCULATING TEMPERATURE DISTRIBUTIONS IN ONE-DIMENSIONAL FREEZING—MELTING AND TEMPERATURE VARIABLE PROPERTY PROBLEMS A.1 General finite difference method. The most flexible method of solving the heat conduction equation for composite bodies with temperature variable pro- perties and general boundary conditions is the method of finite differences. The finite difference equations applying to a general interior node, heat flux boundary condition at x=0 and x=E, and temperature boundary condition at x=O and x=E are presented. A.l.l General interior node. An energy balance for node n, shown in Figure A.l, be- tween two dissimilar materials can be written m+1 m+1 m+1 _ A(n)Tn_l + B(n)Tn + c(n)'rn+1 - D(n) [4.1.1] where the coefficients A(n), B(n), and C(n) are defined by 2kn_i(l-CNBD) A(n) = - 2 (Ax_ + Ax_Ax+) [A.l.l(a)] 112 113 2k l-CNBD 2k l-CN B(n) = 2%? + n+i( ; + g‘&( ED) (Ax_Ax+ +Ax+) (Ax_ +Ax_Ax+) [4.1.100] 2k (l—CNBD) C(n) = - “I 2 [A.l.l(0)] mam. + Aug.) 2kn+&(CNBD) m D(n) = 2 Tn+1 mam. + Ax.) 2k (CNBD) 2k (CNBD + 2%? _ n+i 2 _ 2n-I ) T: (£§x_[§x+ +£§x+) ([5x_ + [5x_[§x+) 2k (CNBD) "‘3 Tm 4.1.1 (1 4' [( AXE + Ax-AX+)] n-l [ ( )] T: is the temperature of node n at time tm, kn++ is the thermal conductivity evaluated between nodes n and n+1, Z§x+ and [3x_ are defined on Figure A41. CNBD is a finite difference parameter that may take on values between zero a and one. If CNBD is set equal to one, the explicit forward difference method results; CNBD being set equal to zero results in the implicit backward difference method. If CNBD is set equal to I the Crank-Nicolson implicit method is obtained. P is the average heat capacity of node n and is given by 114 Interface beineenmateriols a and b material material b AX- AX... O 0 0 n-I n n+l Figure A.l Node geometry and nomenclature for the finite difference solution to the heat conduction equation at the interface between dissimilar materials. AXI II-X)-'~———— >. -——-i-O X 0 q l a 2 Figure A.2 Node geometry and nomenclature for the finite difference solution to the heat conduction equation at the heated surface with a known heat flux. 0 - fixed node X - floating node r—AX'-—.1Ax'+thsz_—sz+—.i l-l l L23 bu I+2 fusion front location Figure A.3 Node geometry and nomenclature for the floating node used in the finite difference solution of the freezing-melting problem. 115 ( I + £3 ( ) P = Ax- P21: + A: POL: [A'l'l(°)] (Pop)_ and (Pcp)+ are the heat capacities corresponding to materials a and b in Figure A.l respectively. [5t is the time step i.e., the difference between times tm and tm+l° If node n is not between two different materials then P simply becomes (f)cp)n, the heat capacity of node n. If the node spacing to the right and left of node n is equal then [3x_ = [§x+. Thus this single equation applies to any interior node, that is, any node not on a boundary. A.l.2 Heat flux boundary condition at x=O. An approximation used by Beck(5) to improve the accuracy of finite difference calculations involving heat flux boundary conditions will be presented here. The method will be refered to as the quarter point method. Referring to Figure A.2, the temperature at point a can be approximated as T: = M? + (1- MT; [4.1.2] The temperature at point a will be used in the approxi- mation for the time derivative instead of the temperature at point 1. An energy balance on node 1 then gives m+1 m+1 _ B(ml + c(1)cr2 _ D(1) [1.1.3] The coefficients B(l), C(1), and D(l) are defined by 116 2(l-CNBD)k]fi "EKL‘ [4.1.300] B(l) (l-AH ) (1- CNBD)kl 0(1) = AtPcE—a - 2 ii [A. 1. 3(b)] [5112 M c )1 (“3ka m 11(1) = - 2 T l: _%?L A"12 ] l + [II-AIIPcPIg + 2(CNBD)kl:}—]Tm At Axl2 1 2(1-CNBD)qm+1 2(CNBD)qm ... Ax]. + All [Aele3(c)] [5x1 is the distance between node 1 and node 2, kl+i is the thermal conductivity evaluated between nodes 1 and 2, and qm and qm+1 are the heat fluxes at times tm and tm+1 r spectively. The parameter A was found by experiment to have a best value of 0.75, thus the term "quarter point method". If q is set to zero, the heat flux boundary con- dition becomes an insulation boundary condition at x=O. A.1.3 Heat flux boundary condition at x=E. The quarter point method is also used at x=E. The re- sulting equation for node N at x=E is 117 Mumfifi + B(N)T§+l = D(N) [1.1.4] The coefficients A(N), B(N), and D(N) are defined by (1- AM pep)N 2(1--0N131))kN_i A(N) = A.1.4( ) At AXNZ I: a] M c ) 2(l-CNBD)k B(N) =.- N + ————15=i 4.1.406) —§—2—. A4- [ I MN) = £82911: _ W Tm At AxN2 N - 2 CNBD k + [(1 ngi‘lm’i + —-—-2—i( ) N' ] TEL]. AXN 2qm+l(l-CNBD) 2qm(CNBD) - [SIN - [SIN [A.1.4(c)] A is again 0.75 and if q=O, we have an insulated boundary at x=E. A.1.4 Temperature boundary condition. For the case in which the boundary temperature is a known function of time we need not write an energy balance for nodes 1 or N, since the new temperature at node 1 or node N is simply equated to the known function. if” = 11(1) [A.l.5(a)] 118 or T§+l = D(N) [A.l.5(b)] where D(1) and D(N) are the known functions evaluated at time tm+l° A.l.5 Solution of the finite difference equations. The result of the energy balances at all N nodes can be written in matrix notation as follows: )- 1? +11 I- 1 3(1) 0(1) 0 0 ....... 0 a? D(1) 1(2) 3(2) 0(2) 0 ....... o 43*1 D(z) '9 4(3) 13(3) cm .....9 8;” = D(3) 0 . ..... . A(N-l) B(N-l) C(N-l) T§:i D(N-l) 0 ............. ... A(N) B(N) T§*1 D(N) _ ..IL. .. I. .- [A.l.6(a)] or in short form m+1 __g = p [A.l.6(b)] The solution for the temperature distribution at time tm+1 is very efficiently obtained by Guass's elimination method. For a discussion of this method see Smith<37). 119 A.2 The solution of freezing-melting problems. A.2.1 Mathematical statement of the problem. The temperature distribution in a substance that melts or freezes isothermally such as a eutectic alloy or water can be described mathematically by(27). 9T 3 33? Pope 3??- = “(ks—87g) [A.2.l] in the solid region, and by 3T 3 ai,’ Pep; 575'?- = 3745. Ex ) [11.2.2] in the liquid region if natural convection can be neglected. At the interface between solid and liquid where E is the position of the interface and T1. is the fusion temperature. An energy balance at the interface be- tween solid and liquid yields1 d€___1|_.___k 9T8 dt’pLsax 3T x= e__" 1‘1 3x 1: 6+] [4.2.4] where fig is the velocity of the interface and pL is the heat of fusion. Equation [A.2.4] couples equations [A.2.1] and [A.2.2] in a highly nonlinear manner. 1. This equation is written for a freezing problem; the interface is moving in the positive x-direction with solid to the left of x=€ and liquid to the right of x=E . Melting \ is a constant deter- mined by solving for the root of equation [2.4.4]. Equation [2.4.4] can be solved easily by the Newton-Raphson iterative method. The time step should be small enough so that El from equation [A.2.22] is much less than x1, the distance from the surface to the first node. The temperature of the boundary node is determined from the given boundary temperature and the temperature of the floating node is fixed at the fusion temperature. The calculation procedure for a single time step is then: 1. Using the fusion front location at time tn, cal- culate the temperature distribution for time tm+1. The body is treated as two separate bodies, one to the left of the fusion front and one to the right; the finite difference equations of section A.l are applied to both. The rear boundary condition of the body to the left of the fusion front is that of constant temperature, i.e. the fusion temperature. The front boundary of the body to the right of the fusion front is also held at the fusion temperature. 2. Calculate the fusion front location for time t m+1 using equations [A.2.l9], [A.2.20], and [A.2.21]. 130 3. Calculate Axl_, Axl+, Ax2_, andAx2+ for use in calculation of the next time step. 4. Return to step 1. Calculation proceeds in this manner until the fusion front crosses the entire body, then the calculation procedure returns to the standard finite difference methods of section A.l. Several test cases were run to check the accuracy of the method and the computer program. Comparison was made with the results of an exact solution for a semi-infinite body initially at a uniform temperature above the fusion temperature with a step change in surface temperature and an approximate solution for a finite slab initially at the fusion temperature with a constant heat flux boundary con- dition at x=O and an insulation boundary condition at x=E. The exact solution is given in Chapter II by equations [2.4.1] through [2.4.9]. A comparison of results from the finite difference solution for both Crank-Nicolson (CNBD=O.5) and backward difference (CNBD=0.0) methods and the exact solution is shown in Figure A.4. The backward difference method seems to be superior to the Crank-Nicolson method for this problem. The approximate solution used to check the accuracy of the computer program was found using an integral energy method similar to that used in Chapter II. The body is considered to be initially completely solid at the fusion temperature; thus, no heat passes the interface between 131 ICC 0°75 0°50 +.. LT 4352 —T +.. TT 050 025 - exact solution 1 1 O backward difference approximatlon T 0 Crank-Nicolson approximation j‘T DIMENSIONLESS TEMPERATURE, T = (T(x,1)'T(O,t))/(T(x,0)-“0,1” Y 0.00 J J J J J__1_1 J J _l __l __l J J l J L J _JJ_J 00 I6 so 46 so omeusrowLess TIME, Ting/x2 Figure A.4 A comparison of finite difference solutions to an exact solution of the freezing-melting problem. 132 solid and liquid and only the prOperties of the liquid and the heat of fusion become involved in the solution. After the fusion front moves across the entire body, an exact solution for a finite slab with a constant heat flux at x=O and insulated at x=E was used to describe the tem- perature distribution. The solution is T-T 2 T2 " un7—kx' ‘ A[E E] " B[E E] ' 757w O‘E‘E [A.2.23(a)] - T'Tf , ‘ €‘x‘ Ts ... M; = o, 1-7m, 13-3- 1 [A.2.23(b)] where A and B are given by L+ L+ L+2 % A = .. ___Q_+ -—-£— 11.2.23 —% [(6/13) 4(6/E)2] [ (0)] and L+ L+ L+2 % z ,1, l + ___9_]- .3 __9._ _I.‘_c:__ (e/E) 2(e/E) (es/N 4(e/E)4 [A.2.23(d)] The location of the fusion front e/E is found by solving the algebraic equation as Mm we as) warm m] The time ‘Th is the time for the entire body to melt and is given by arm = EL; +% [1,; + %L;2T + "16'- [A.2.25] 133 After the entire body melts (’T>Tm), the solution is given by T-T 2 "' .. f .. a _ 1‘. _ Tz—mi-Tt-fhz) E A+B+% co 2 2 2 2B - l -n 1T 1' _ g + _7T2 ———_n2 e cos[n1r(l En n:l [A.2.26] Results of the approximate solution and the finite difference solution using the backward difference ap— proximation are shown in Figure A.5. There seems to be quite good agreement between the two solutions. 134 I l-500 finife difference r- — infeqrai oppro’ximafion I425 r- L’,‘ -2-o, Taco-o 0750 0375 DIMENSIONLESS TEMPERATURE, (T-Tfi/(qE/l‘) 0°00 0° 075 omensnowuzss TIME, range/E: Figure A.5 A comparison of the finite difference solution to an approximate solution of the freezing- melting problem. APPENDIX B THE PARAMETER ESTIMATION PROCEDURE AND THE CALCULATION OF SENSITIVITY COEFFICIENTS 3.1 The Gauss-Newton or linearization method of minimizing the sum of squares function The sum of squares function 18(40) P(E) = [z - w<2>]y‘1[z - my] [13.1.1] Where 1 is the vector of observations of n thermocouples and m time steps. r- '1 (m x 1) ll Y(nm x l) Y = .prz) [B.l.3(a)] If we substitute equation [8.1.2] into equation [3.1.3], we obtain £3?le - 1(2°)] + z(p°)<;3_* - p°> = 9 [3.1.4] Solving for 2* in equation [B.l.4] 3* = 13° + N’l;'(§°)_Y_‘l [x - 3(g°)] [13.1.5] 1. Note the prime (') denotes transposition. 138 where the square matrix K is defined by 1; = g'l’l; [B.l.6] The true value of 2* is usually not reached on the first iteration so the iterative formula 2" = 2“ + i‘lz'<2““1>1‘1[1 - mk'H] [3.1.7] is applied. The iteration can be stopped when all the pi - pE-l have satisfied some predetermined criterion such as k k-l p < --E--' 0.001 B.2 Calculation of sensitivity coefficients The matrix Q is called the matrix of sensitivity coefficients. The sensitivity coefficients are an important part of parameter estimation, not only in the Gauss-Newton algorithm but in the determination of Optimum experiments. The numerical approximation of these derivatives is not difficult once the mathematical model has been programmed; the differentiation can be approximated as 3W(£:Xioti) EDPJ w(plo'°°9pl+ 8p19°npk9xi9ti) " w(p1!°-01p1v4tgovxivt1) 8p [13.2.1] J' 139 where 8 is some small value, say 0.001. In order to test the method and the computer program, sensitivity coefficients for a test case were computed using the computer program and compared with sensitivity coefficients calculated using the exact solution. The test case used a finite body exposed to a constant heat flux at one surface (x=0) and insulated at the other (x=E). The exact solution for the dimensionless sensitivity coefficients is 0° .7.“ 32 = .. LL—E-x 2 - 32 + —a— .2— 9-..?" 27564;”) qE k Bk: 6E2 '1r2~n= n2 E 09 —n2'n'21' + 21’ e cos(n1T x) [13.2.2] n=l E a. 2 2 qE k —(--———73 33c a T[l + ZZ e-n ’fl' Tcos (n1ETx] [R25] p n=1 The results are shown in Figure B.l and Figure B.2; the numerical approximation to the derivatives seems to work quite well for the case of no change of phase. 140 042» I-O 04.. OVTS // f/ ’ 0'5 o-o ’2 "/ -04 >- 0.25 -- exact solution _, 0,2 __ o o finite difference solution .. 0.3 ._ .. 0.4 l J 1 1 J 04) 04 ()2 (>3 (>4 (15 'r sat/£2 Figure B.l A comparison of exact and numerical ap- proximation of the dimensionless sensitivity coefficient SE(x,t) for a finite body with constant heat flux, q, at x=0 and insulated at x=E. ..(y5p -O~4- -0" .— 141 — exact solution e e finite difference solution 0.0 , I l L j o-o O-l (>2 (>3 O~4 (>5 7 I at IE2 Figure B.2 A comparison of exact and numerical ap- proximation of the dimensionless sensitivity coefficient 35° (x,t) for a finite body with constant heat flux, q, at x=0 and insulated at x=E. APPENDIX C THE CALCULATION OF SURFACE HEAT FLUX FROM AN. INTERNAL TEMPERATURE HISTORY The problem of finding the surface heat flux given the temperature history of a body at some internal point is known as the inverse heat conduction problem. That is, the boundary condition at one surface (the heat flux) is un- known, while the temperature (or solution to the normal boundary value problem) at some internal point is known. The heat conduction equation for a one dimensional homogeneous body is {2(1‘33— k-¥=) pcpi- at [0.1] and the boundary and initial conditions for the inverse problem, assuming for simplicity the body to be insulated at x=E, and the temperature history for some point xse is known, are T(e,t) . T.(t) [0.2] 1'33: a 0 V [0.3] T(x,o) = 11(1) [c.4] The problem is to calculate q(0,t) . -1: Jig-2&3)- [c.5] 142 143 For a review of the various methods, the reader is referred to the paper by Beck<35). To illustrate the diffi- culty encountered by most methods of solving the inverse problem, Beck considered the effect of a short heat pulse of duration ‘Tp at x=0 upon the temperature rise at x=E, the insulated surface of the body. The ratio of the tem- perature rise to the maximum temperature rise [T(E,t)-T1]/ [T(E)max'Ti] versus T/‘T p for several values of 1p is shown in Figure C.l. Note that the ratio of temperature rise to maximum temperature rise at time ‘T/“Tp=l becomes very small as ‘Tp decreases below 0.1. Thus very little information about the heat pulse reaches the insulated surface in time 7p for 1p less than 0.1. As T/‘Tp becomes larger than 1; however, more information becomes available about the heat pulse. Beck used this fact in his concept of "future temperatures". That is, to calculate the average heat flux over the time interval ‘Tp he employed temperatures at times 2‘Tp, 3‘Tp,... The temperatures at times 2‘Tp, 31'p,... are called future temperatures. In general, the smaller‘l’p the more future temperatures are needed. The procedure to solve the inverse problem for an arbitrarily varying heat flux is to divide the heat flux into a finite number of heat pulses. The duration of pulse qm is from time tm-l to time tm. This is illustrated in Figure 0.2. The duration of the various pulses need not be same. ql is found by minimizing the sum of squares function ~l .~ . .... . v. Vg'vrw' ’ 144 '1': (>8? . (325 04 0'5 " o-os o-4 - o-oz Tifil- TL “(El - TI tn, Figure C.1 Ratio of temperature rise to maximum temperature rise versus ratio of dimen- sionless time to dimensionless duration of heat pulse at the insulated surface of a plate with a heat pulse applied to the other surface. 0 HEAT FLUX 4D .0 03 N ‘ i i I , I Iez a 3 f I I I I T I I '0 'I ‘2 '3 'm-I 'm TIME Figure 0.2 Illustration of how a time varying heat flux is broken up into a series of heat pulses. 145 I 2 _ l+i l+i F(qp - 2119 8L g 0 .I(>6~ u. E 0'4? I32 I (yak I 03:0” 6:2" 034' 0-6 as m l-2 l-4 omensnomess nuance Figure 0.4 Heat flux calculated from noisy data for 2 future temperatures for a slab with a step in flux at one surface and insulated at the other surface. lllilif 93 03177 6226 12 3 I" L H v. ."I H H " III l H H H ”I H II