THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE OPTIMIZATION OF DOMESTIC HEATING SYSTEMS Adel H, Eltims ahy..'.:.' 0 ~ A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan 1967 April, 1967 IP-779

Doctoral Committee: Professor Louis F. Kazda, Chairman Associate Professor Keki B. Irani Dr. Vernon Lo Larrowe, Research Engineer Associate Professor Herman Merte, Jr. Associate Professor William A. Porter

ACKNOWLEDGMENTS The author wishes to express his sincere appreciation to the chairman of his doctoral committee, Professor L.F. Kazda, for his advice and counsel during this study. He is also grateful to the other members of the committee, Professors K.B. Irani, V.L. Larrowe, W.A. Porter, and especially H. Merte for their comments and suggestions. He also wishes to express his thanks to the staff of the Maxitrol Company for offering their technical assistance and valuable information in obtaining data for this study. The author is indebted to the Industry Program of the College of Engineering for the reproduction of the dissertation. iii

TABLE OF CONTENTS Page ACKNOWI0DGMENTS.............. iii LIST OF TABLES............................... vii LIST OF FIGURES.................................. viii LIST OF APPENDICES............................................ xiv I INTRODUCTION.......................................... 1 1.1 Outline of the Problem.......... 1 1.2 Work Relating to the Problem.................. 4 1.3 Comfort........ 6 1.4 Components of the Domestic Heating System...... 8 II THE MATHEMATICAL MODEL OF THE ROOM...........o.... 12 2.1 The Physical Model................ 12 2.2 The Room Spac 13 2.2 The Room Space................................... 13 2.2.1 Air and Heat Distribution in the Room Space...,..................... 15 2.2.2 Region 1 or the Jet Region Model........ 17 2.2.3 Region 2 Model......................... 27 2.2.4 Region 3 Model................... 5 30 2.3 The Wall Model................................ 32 2. 4 The Room System Model......................... 39 III THE EXPERIMENTAL CHECK OF THE ROOM MODEL............ 43 3.1 System Investigated........................... 43 3.1.1 General Description..................... 43 3.1.2 Construction Details.................... 46 3.1-3 Instrumentation......................... 46 3.1.3.1 Temperature Measurements....... 46 3.1.3.2 Velocity and Rate of Flow Measurements..................... 49 3.2 Calculation of System Parameters.............. 49 3.2.1 Specific Heat of Air................... 49 3.2.2 Density of Air.......................... 49 iv

TABLE OF CONTENTS (CONT'D) Page 3.2.3 Volume Calculations................... 49 3.2.4 Thermal Conductivities.................. 51 3.2.5 Thermal Resistances.................... 51 3.2.6 Equivalent Capacity of the Room Boundaries 51 3.3 The Experimental Test................... 55 3.4 Computer Simulation. 58 3.5 Evaluation of Model........................... 64 IV THE CONVENTIONAL HEATING SYSTEM ANALYSIS................ 72 4.1 The Mathematical Model of the Controller...... 72 4.2 The Mathematical Model of the Heat Exchanger... 74 4.3 The Mathematical Model of the Hot Air Duct..... 81 4.4 The Mathematical Model of the Gas Control Valve 83 4.5 Analog Computer Model.......................... 83 4.6 The Heating System Analysis.................. 86 4.6.1 Effect of Thermostat Hysteresis on Cycling Amplitude and Cycling Period.... 88 4.6.2 Effect of Thermostat Time Constant on Cycling Amplitude and Period............ 92 4.6.3 Effect of Furnace Time Constant on Cycling Amplitude and Cycling Period.... 92 4.7 Evaluation of Results...9.................... 95 V OPTIMIZATION OF THE HEATING SYSTEM................... 99 5.1 Reformulation of the Mathematical Model........ 99 5.2 Formulation of the Variational System......... 102 5.3 The Optimization Criterion......... 105 5.4 The Optimal Control Law.................... 106 VI THE ANALOG COMPUTER SIMULATION OF THE OPTIMAL HEATING SYSTEM........ o....... o.............................. 118 6.1 The System Equations......................... 118 6.2'The Parameters k31, k32, and k33 of the Optimum Controller............................. 120 6.3 Analog Computer Simulation of the Optimum System.............o.o.......................... 121 6.3.1 Fixed Part of the Heating System..,..... 121 6.3.2 The Simulation of the Optimum Controller 123 6.3.3 Equilibrium Values..................... 123 v

TABLE OF CONTENTS (CONT'D) Page 6.4 Analysis of the Optimum Heating System.......... 124 6.5 Conclusion...................................... 131 VII A SUB-OPTIMAL SYSTEM...................................... 135 7.1 The Proposed Sub-Optimum Heating System......... 136 7.2 Analysis of the Sub-Optimum Heating System...... 140 7.3 The Practical Implementation of the Sub-Optimal Controller..................................... 148 Summary............................ 153 VIII CONCLUSION................................. 154 REFERENCES............................................ 183 vi

LIST OF TABLES Table Page 3.I Computer Potentiometer Settings of Figure 3.16..... 64 C-I Properties of the Building Materials................ 168 D-I a. Wall Properties............................ 174 b. Ceiling Properties.................. 175 c. Floor Properties.................... 175 vii

LIST OF FIGURES Figure Page 1.1 Components of the Heating System.................11 2.1 Side View of Room Showing the Different Regions..... 14 2.2 Plan View of the Room............................... 14 2.3 Heat Transfer Near the Wall.................. 16 2.4 Side View of the Room Showing Heat Transfer,..... 18 2.5 Side View of the Room Space Showing Air Distribution 18 2.6 The Three Phases of the Jet...................... 20 2.7 Velocity Profile for a Free Jet.............. 20 2.8 Velocity Distribution at a Cross-Section Normal to to the Axis of the Jet........................ 24 2.9 The Control Volume in the Jet Region Corresponding to a Distance x.................................... 24 2.10 Control Volume of Region 1....................... 24 2.11 Variation of Tx with x........................ 26 2.12 State Variable Diagram for the Jet Region.......... 28 2.13 Control Volume of Region 2.......................... 28 2.14 State Variable Diagram for the Second Region....n... 31 2.15 Control Volume of Region 3*........................ 31 2.16 State Variable Diagram for the Third Region......... 31 2.17 The Wall Thermal Circuit............................ 35 2.18 State Variable Diagram for Part of the Wall Facing Region 2...a... a................................ 40 2.19 State Variable Diagram for Part of the Wall Facing Region 3............................ 40 2.20 State Variable Diagram for the Room System.......... 42 viii

LIST OF FIGURES (CONT'D) Figure Page 3.1 View of Interior of the Maxitrol Room.................... 44 3.2 The Maxitrol Test Room with Dimenstion Shown........ 44 3.3 Construction Details of-the Test Room Boundaries..... 47 3.4 The Setup for Temperature Measurement.............. 48 3.5 The Test Room with Regions shown.....................50 3.6 Input Temperature for Capacity Measurement.5..... 53.. 3.7 Average Temperature of Region 3 for Capacity Measurement.....................,................. 53 3.8 Inside Surface Wall Temperature for Capacity Measurement.............................. 53 3.9 Outside Surface Wall Temperature for Capacity Measurement.................5.............................. 53 3.10 Location of Points at which Temperatures were taken during the Test..................................... 57 3.11 Location of Points on the Inside Surface of the South Wall, where Temperatures were Recorded During the Test.................................................. 60 3.12 Piecewise Linearization of the Input'Function........ 61 3.13 The Input Temperature................................. 62 3.14 Analog Computer Diagram to Generate the Function Ti... A; *.e * * *.e ~ --.,,,, ~~....... e'e-.~~.~.. ~~ ~ ~....~,' ~~~~ ~ ~ 3.15 Analog Computer Diagram for the Room System........... 65 3.16 Average Temperature in Region 1...................... 66 3.17 Average Temperature in Region 2.................... 67 3.18 Average Temperature in Region 3............67 3.19 Average Wall Surface Temperature Facing Region 2...... 68 3.20 Average Wall Surface Temperature Facing Region 3..... 68 ix

LIST OF FIGURES (CONT'D) Figure Page 4.1 A Bimetallic Thermostat.............................. 75 4.2 State Variable Diagram of the Thermostat............. 75 4.3 Hysteresis Function.................................. 75 4.4 A Simple Diagram for a Typical Gas Furnace........... 77 4.5 Heat Transfer in the Furnace......................... 78 4.6 State Variable Diagram for the Heat Exchanger........ 80 4.7 Schematic of an Air Duct............................. 82 4.8 State Variable Diagram for the Hot Air Duct.......... 84 4.9 Analog Computer Diagram for the Thermostat........... 85 4.10 Analog Computer Diagram for the Heat Exchanger....... 85 4.11 Analog Computer Simulation for the Simplified Room Model............................................... 87 4.12 A Sample of Recordings.............................. 89 4.13 A Sample of Recordings............................... 89 4.14 Conventional System's Cycling Amplitude-Hysteresis Curves with TT as a Parameter...................... 90 4.15 Conventional System's Cycling-Period-Hysteresis Curves with TT as a Parameter...................... 90 4.16 Conventional System's Cycling Amplitude-Hysteresis Curves with TF as a Parameter..................... 91 4.17 Conventional System's Cycling Period-Hysteresis Curves with Tp as a Parameter..................... 91 4.18 Conventional System's Cycling Amplitude -TT Curves with q as a Parameter.............................. 93 4.19 Conventional System's Cycling Period -TT Curves with q as a Parameter................................... 93 4.20 A Sample of Recordings............................... 94 x

LIST OF FIGURES (CONT'D) Figure Page 4.21 A Sample of Recordings.............................. 94 4.22 A Sample of Recordings,................. 94 4.23 Conventional System's Cycling Amplitude -TF Curves with q as a Parameter............................ 96 4.24 Conventional System's Cycling Period -TF Curves with q as a Parameter........................... 96 5.1 Block Diagram of the Optimum Heating System.......... 117 6.1 Block Diagram of the Optimum Heating System......... 122 6.2 Analog Computer Simulation for the Generation of the Equilibrium Values........................... 125 6.3 Analog Computer Diagram for the Optimum Heating System 126 6.4 Room Temperature Variation Response for Case 1...... 128 6.5 Surface Wall Temperature Variation Response for Case 1..................................... 128 6.6 Heat Exchanger Temperature Variation Response for Case 1............................ 128 6.7 Flame Temperature Variation Response for Case 1...... 128 6.8 Room Temperature Variation Response for Case 2....... 130 6.9 Surface Wall Temperature Variation Response for Case 2......................... 130 6.10 Heat Exchanger Temperature Variation Response for Case 2............................................... 130 6.11 Flame Temperature Variation Response for Case 2...... 130 6.12 Room Temperature Variation Response for Case 3....... 132 6.13 Surface Wall Temperature Variation Response for Case 3............................................... 132 6.14 Heat Exchanger Temperature Variation Response for Case 3.................................... 132 xi

LIST OF FIGURES (CONT'D) Figure Page 6.15 Flame Temperature Variation Response for Case 3.... 132 7.1 Block Diagram Schematic of the Sub-Optimum Heating System.................................... 139 7.2 Cycling Amplitude-kw Plots with q as a Parameter... 141 7.3 Cycling Period-kw Plots with q as a Parameter...... 141 7.4 Cycling Amplitude-kw Plots with TF as a Parameter.. 143 7.5 Cycling Period-kw Plots with TF as a Parameter..... 143 7.6 Cycling Amplitude-kw Plots with TT as a Parameter.. 144 7.7 Cycling Period-kw Plots with TT as a Parameter..... 144 7.8 A Sample of Recordings............................ 145 7.9 A Sample of Recordings............................ 145 7.10 The Solid State Sub-Optimal Controller............. 150 7.11 Block Diagram of ths Sub-Optimal Controller.........151 A-1 Configuration Showing the Different Regions........ 157 B-l An Element of Solid with Heat Flow in one Direction 161 B-2 Transmission Line Lumped Model..................... 163 C-l The Experimental Set-up for Capacity Estimation.... 167 C-2 Construction and Electrical Model of the Wall...... 170 C-3 Construction and Electrical Model of the Floor..... 170 C-4 Construction and Electrical Model for the Ceiling.. 170 D-1 Dimensions of Region 1............................. 173 D-2 Temperatures of Points in Region 1................. 178 D-3 Temperatures of Points in Region 2................ 179 xii

LIST OF FIGURES (CONT'D) Figure Page D-4 Temperatures of Points in Region 3.................. 180 D-5 Surface Temperatures of Points Facing Region 2...... 181 D-6 Surface Temperatures of Points Facing Region 3...... 182 xiii

LIST OF APPENDICES Appendix Page A TIME DELAY CALCULATIONS............................ 157 A-i Time Delay in Region 1......................... 157 A-2 Time Delay in Region 2......................... 157 A-3 Time Delay in Region 3..................... 158 B BASICS OF THERMAL CIRCUITS.......................... 160 B-1 Flow of Heat in Solids......................... 160 B-2 Flow of Charge in a Noninductive Transmission Linein.................................... 163 B-3 The Analogy between Thermal and Electrical Constants and Variables...................... 165 B-4 Scale Problems...................6....... 16 C HEAT CAPACITY ESTIMATION OF-THE ROOM BOUNDARY....... 167 C-l Estimation of the Capacity of Region TwoBoundary......................... 168 C-2 Estimation of the Heat Capacity of Region Three-Boundary........................... 169 D EXPERIMENTAL MODEL CALCULATIONS AND TEMPERATURE PLOTS............................................... 173 D-l Volume Calculations........................ 173 D-2 Thermal Conductivity Calculations............. 173 D-3 Thermal Resistance Calculations............... 176 xiv

CHAPTER I INTRODUCTION 1.1 Outline of the Problem A review of the research and work which has been done, and is being done, in the area of gas-fired forced-air controlled comfort heating systems reveals little documented effort expended to study the heating system from a control system point of view. Furthermore practically no effort has been expended in trying to optimize these systems using the techniques of modern control system theory. By performing a library research and a study of the pertinent literature in the heating systems area, it is found that two approaches to this problem are available. a) One approach is to use full Sacale experiments and utilize a testing procedure designed to indicate the effect of a given system modification. These experiments are practicable in principle, but they are exceedingly difficult to realize except under controlled laboratory conditions. In addition, they are very time consuming and present a very lengthy and expensive program. This approach has been used for many years and has provided many useful results for the warm air heating system designer. It is especially useful when the objective is to investigate design improvements for particular components of the heating system. Examples of this type of endeavor are found in research bulletins issued by the Engineering Experiment Station of the University of Illinois and by the National Warm Air Heating and Air Conditioning Association. -1

-2b) The second approach is to utilize an analytical model to represent the system. The advantages of the analytical approach are relative inexpense, ease of duplication, and the relative ease of evaluating results. Furthermore, the analytical approach often provides valuable insight which leads to more complete understanding of the actual system. It must be emphasized, of course, that the accuracy of the results obtained utilizing the mathematical model is entirely dependent on the degree to which the model represents the actual system. In the past few years other researchers have performed analytical studies of the domestic heating process. However, these studies were primarily limited to open loop, steady state situations,(4,5) and were not concerned with the dynamic response of the system. Incentive was thus provided to study the problem which is defined as follows: 1o To describe the dynamic behavior of the domestic room system using the knowledge of principles of thermodynamics and air flow dynamics. Previous experimental results will also be used where applicable to assist in describing the behavior of the room system. 2. To formulate an appropriate mathematical model using the results in (1) to describe the dynamic operation of the room system. This will be an approximate distributed parameter forced-air room model or a lumped parameter model, (1) and (2) are treated in Chapter IIo 30 To simulate the mathematical model of the room so derived on an electronic differential analyzer and to compare the results obtained with experimental data since a thorough evaluation of the

-3applicability of the model would require a parallel comparison between a real installation and the room model. This experimental check is presented in Chapter III. 4. From a knowledge of the physics of the various elements that compose the rest of a typical heating system, mathematical models are to be determined that would adequately define the dynamic behavior of each of the elements in their pertinent operating ranges. The mathematical models of the elements thus derived together with the room's mathematical model are then to be simulated on an electronic differential analyser. 5. Typical heating system parameter variations are to be studied and analyzed. Results are then to be related to experimental observations in the heating systems industry. Results of analyses are to be used also for the purpose of comparison with other systems. (4) and (5) are treated in Chapter IV. 6. To apply modern control theory methods of analysis to the gas-fired forced-air heating system model with the purpose of providing optimum comfort and cost. The control problem is as follows: Given a fixed plant which consists of a furnace, a duct system, a domestic space and the walls, it is required to find the optimum controlling unit that best controls the heating process according to a specified criterion. This is presented in Chapter V. 7. To simulate the optimum heating system developed on the electronic analog computer for the purpose of comparison with the conventional heating system and other systems. Such analysis is shown in Chapter VI.

-48. To study the optimum controller from the practical point of view and if necessary devise a sub-optimum controller using properties of both the optimum controller and the conventional controller. 9. To simulate the sub-optimum heating system on the electronic analog computer in order to compare it with the optimum and conventional heating systems. 10. Knowing the dynamic characteristics of the sub-optimum controller, a component possessing these characteristics, and the method of implementing it practically is presented. (8), (9), and (10) are presented in Chapter VII. 1.2 Work Relating to the Problem In work currently being carried out by the University of Illinois, sponsored by the American Boiler Manufacturers, a model house is being tested. The main objective of the test is to observe temperature patterns in the house. The walls and air space of the model house are instrumented with 200 thermocouples, 125 of which are used to record temperature manually using a type K-l potentiometer as the measuring instrument. Once the data is recorded, correlation routines are used to determine variances, means, etc. and thus establish temperature patterns. Temperatures are recorded at approximately 90", 60", 30", and 3" intervals from ceiling to floor. Readings are taken at 8 a.m., 11 a.m., 4 p.m., and 10 p.m. every day. The data has recently been put on IBM cards for processing. In another project at the University of Illinois which is sponsored by the National Warm Air Heating Associationwork is being done on proper duct sizing for heating and cooling systems. The objective of the project is to assess the effect of joints, bends, elbows, etc. in order to properly

-5size the ducts and furnaces for residential structures. In addition studies are also being conducted in order to determine the effects of air filters, humidity, air velocity, sound of air passing through ducts, and how the air is introduced to the room. In order to conduct these tests a room 18' x 13' x 15 in size is avialable whose outside temperatures can be controlled. The Maxitrol Company in Detroit is currently doing experimental research work on modulation systems, (i.e. systems in which the heat output of the furnace is controlled in a continuous fashion) and in determining the factors of major importance involved in the design of these systems. It has been found that most domestic gas-fired furnaces are oversized and that design specification should be modified. For this reason a thermostatically controlled enclosure 12' x 12' x 8' in size has been constructed. It is designed to give uniform air temperature distribution throughout the structure. At the University of Wisconsin their studies to date are in the industrial systems area. Their research work has been sponsored by ASHRAE. The objective of their work is to introduce automatic control concepts into the ASHRAE literature in order to familiarize the membership with these more theoretical concepts. 5', They are also studying the effects of thermostats on the control of comfort in home heating systems using a simulator. Provision was made to vary thermostat delay, thermostat differential, and outdoor temperature. They are also studying the heating anticipation problem for their out-of-door-in-door thermostat system of control. At the University of California in Los Angeles, the thermal behavior of a one room house subject to daily variations was studied. This

-6work is the result of research sponsored by the American Society of Heating and Air Conditioning Engineers Inc..(9 The thermal-electrical analogy was used. Input devices were developed to represent diurnal variations in air temperature, solar energy absorbed, and net long wave radiation exchange. In parallel with this study, research work on the optimal control of gas-fired warm air heating systems sponsored by the American Gas Association is being conducted at the Electrical Engineering Department, University (12) of Michigan. (2) The main objective of the research is to apply modern control theory methods of analysis to gas-fired home comfort-producing equipment with the purpose of providing superior comfort, and better system integration of components. The objective also is to define one or more typical gas-fired, controlled comfort systems, which can be used as guidelines in the analysis and design of most residential heating and or cooling systems. Finally to simulate typical design conditions on an electronic differential analyser in order to determine the optimum design parameters. 1.3 Comfort The main objective of the gas-fired forced air domestic heating industry is to produce a heating system which provides. the best human comfort for the consumer while operating within such constraints as economy, safety, and reliability. It is evident that any attempts towards a realization of this goal requires a consideration of many factors, one of the most important being the determination and definition of the highly subjective quantity:.known as human comfort....

As difficult as this problem is, however, some reasonable conclusions may be reached, with regard to human comfort, by investigating those factors which affect the heat balance of the human body. Heat exchange between the human body and its environment occurs in three ways: conduction, convection, and radiation. For this reason the most important environmental factors affecting the human body heat balance relation consist of: air temperature, floor surface temperature, air-water content, air velocity in the vicinity of the human body, and temperatures of material objects in contact with the human body. The latter factor, namely contact heat loss, can generally be neglected since uncomfortable contact is normally avoided. Investigations(6) have indicated that the following relationships generally exist between the major environmental factors and human comfort. 1. Air Temperature Air temperature is the most important and most commonly recognized factor which contributes to comfort, An acceptable range is 700 to 77~F at any point within the occupied zone and at any time.(6) 2. Floor Surface Temperatures Floor surface temperatures which are either too high or too low have a great effect on comfort. A heating system may be considered to have excellent performance with regard to floor surface temperatures if temperatures no lower than 650F are maintained on the exposed areas of the floor surface. Furthermore floor surface temperatures should not exceed approximately 80~F at any time.

-83. Air Motion Air motion is an important factor affecting human comfort. Although small values of air velocity will be tolerated, large values will not be because of the general physical discomfort and excessive convection (7) heat loss.7 4. Humidity For health considerations high values of relative humidity are (7) not tolerated. It should not exceed 60 percent, at any part of the occupied zone. 5. Filtered Air For health considerations it is necessary to keep the air clean at all times. In practice the motion of the air is controlled by the speed of the blower, the humidity in the system limited by using a humidifier, and the air is filtered by the addition of an air-filter to the system. With this in mind, it is evident that the room temperature has perhaps the greatest influence on human comfort, and consequently is considered to be the primary criterion for judging systems improvement. Of course throughout the study other factors such as systems cost, and power limitations will also be considered. 1.4 Components of the Domestic Heating System In order to study the control problems associated with a forcedair gas-fired domestic heating system it is necessary to investigate the various elements which constitute the present heating system. Any suitable listing of these elements must include the following basic group:

-91. Domestic Space The region to be comfort controlled. 2. Room Boundaries It includes walls, ceiling, and floor. 3. Furnace The system component which supplies the thermal energy. 4. Controller (Room Thermostat) The logic unit which supplies signals to the furnace to regulate the amount of thermal energy supplied to the domestic space. 5. Hot Air Duct The system component which transfers the heated air from the furnace to the domestic space. 6. Cold Air Duct. The system component which transfers the air from the domestic space to the furnace. 7. Supply and Return Air Outlets The system components through which the air is introduced into and out of the domestic space respectively. 8. Blower The component which circulates the air throughout the system. 9. Humidifier The component which regulates the percentage of moisture in the circulating air. 10. Gas Control Valve The system component which releases combustable gas to the furnace heat exchanger on signal from the controller unit.

-1011. Air Filter The component which removes lint and dust particlas from the air circulated throughout the heating system. In general all of the above elements affect the dynamic perfromance of the system; effects of some of the components as will be seen, however, are negligable. It should also be noted at this time that in the present day heating systems the controller (listed above) generally consists of a simple bimetallic element thermostat. As indicated, however, one of the objectives of this thesis is to determine if other possible and better controllers exist. The above listed components may be integrated into a schematic block diagram indicating their relationship to the overall system. Such a diagram is shown in Figure 1.1.

Outside Temperature oom Boundaries Reference Gas Contro ot Air ^ >ThermostatH Valv Furnace Blower HumidifierH Iui Room Space~Temperature e G __ IAir Cold Air Falter I Duct Figure 1.1 Components of the Heating System.

CHAPTER II THE MATHEMATICAL MODEL OF THE ROOM The microscopic room model as discussed herein refers to the study of the physics of a domestic room space under the conditions of thermal forced air convection in which the room space is not assumed to be at one uniform temperature. Realizing that the room is a complex distributed parameter system, the object of this chapter is to develop a lumped parameter mathematical model which possesses the predominate features of the distributed parameter system but which is amendable to computer simulation studies. In the following section a physical model for the room is assumed which will be used to derive the system mathematical model. Then the room is divided into two main parts: the room space and the walls, ceiling and floor. From a knowledge of the physics of the various elements that compose the room system, the mathematical model of each element is determined and then combined together to give the total model. The total mathematical model is developed such that it would adequately define the dynamic behavior in the pertinent operating range. 2.1 The Physical Model The first problem that arises in this mathematical modeling is to choose a physical model and then analyze the model using the principles of thermodynamics and air flow dynamics. Let us consider a room with height H, width W, and length L. Let the room have three inside walls and one outside wall. An outside wall is a wall which possesses an exterior exposure to the climatic -12

-13 - elements. Let the rest of the house be at a temperature range which is equal to the temperature range of this room. This means that there will be no heat transfer through the three inside walls or the floor or the ceiling. Heat will transfer to the outside only through the outside wall. Let the air inlet to the room be a rectangular opening and let the jet discharge parallel to the inside wall with one edge of the outlet coinciding with the inside wall. See Figures 2.1 and 2.2. To simplify the analysis of the system, no windows or doors are assumed. It will be seen, however, in the next chapters that this does not change the concepts obtained from analysis, since this is essentially equivalent to changing the system parameters; conductance and capacitance. 2.2 The Room Space The usual approach in the study of distributed parameter systems is to replace the continuous space variables by a set of discrete space variables. The accuracy required determines the number of discrete steps needed. For high accuracy the number of steps will be small, while for lower accuracy fewer steps are required. In the approach developed herein uses thermodynamic principles and the air flow pattern in a confined space. The space is divided into regions each having (a) a uniform temperature over-all or (b) having a temperature which varies only in one dimension. In the analysis, one of the regions will be used as a reference region and is also used to control the comfort of space. This will be called the third region.

-14Ceiling 2 I I \ I ~c3 l ^.l A ^~~~~~~~~~~~~~~~~~~~! _.!:,id e r O i et, ins i\e wal. rd W Fu 2 Pa i oh o 0 F.lo or Figure 2,1. Side View of Room Showing th Different Regions' inside wailFigure 2.2. Plan View of the Room.

-152.2.1 Air and Heat Distribution in the Room Space Before going into the details of mathematical modeling, a qualitative description of the air and heat distribution in the room is necessary. The heat transfer involved with the room is by thermal conduction, thermal convection, or thermal radiation. Heat transfer by thermal conduction results due to having molecules of high kinetic energy and molecules of lower kinetic energy. The molecules of the high kinetic energy will transmit part of their energy to the ones with lower kinetic energy. But since the temperature is proportional to the average kinetic energy, thermal transfer will be in the direction of decreasing temperature. In the wall the significant mechanism of heat transfer will be always thermal conduction. If, in addition to conduction, energy transfer occurs by eddy mixing and diffusion it is called thermal convection. In Figure 2.3, in the laminar layer adjacent to the wall, heat transfer occurs by conduction. In the transition region, eddy mixing and conduction effects are significant. In the eddy region heat transfer occurs mainly by eddy mixing. When the air currents are produced by sources external to the heat transfer region, it is called forced convection. If the air currents are generated internally, as a result of non-homogenuous densities arising from the temperature variations, it is called free convection. Both kinds of convection are occurring in the room. In heat transfer by radiation, the internal energy of the source is transferred to electromagnetic energy for transmission, then back to internal energy at the receiver.

tr -Outside 5D Wall R ) kI X >^ Laminar Region F- (C Transition Region Eddy Region (D e ~ o o o.o* (DY.0oD" o/. o0 c+ Room Inside

-17Heat transfer by radiation will be neglected with respect to the quantities transferred by conduction and convection. This is quite justified since no radiating objects are assumed in the room and the quantity of heat radiated from the sun will be included in an outside generator. Figure 2.4 shows how the heat is transferred. Since the jet is adjacent to the wall when it expands into the room, its form will be a semi cone shape. The temperature of an element volume of air coming out from the register will decrease gradually according to a certain variation mainly due to the cooler air already present in the room. At a certain height (a character of the jet and the height of the room which will be specified later), the semi-conical form of the jet disappears and air circulation in the room starts. Based on some experiments done before on the air distribution in the room,(l) the following regions (shown in Figure 2.5) will be used to establish the mathematical model. The choice of the three regions as seen is based on the way the air distributes itself in the room. In the following sections each of the three regions will be studied in detail. 2.2.2 Region 1 or the Jet Region Model In practical room air distribution the airstream is not free since it is influeced by walls, ceilings, floors, and obstruction. How(2) ever, tests have shown that(2) in cases where the outlet area itself is small compared to the dimensions of the space normal to the jet, the jet may be considered free as long as: x < 1.3

-18Ceiling No Heat Transfer O O CH0) o 0 rA IIO cd cd ~r _ O ed 0 No Heat Transfer Floor Figure 2.4. Side View of the Room Showing Heat Transfer. Figure 2.5. Side View of the Room Space Showing Air Distribution.

-19where x = Distance from face of outlet in feet As = Cross sectional area of the confined space in square feet. We shall assume that this inequality is true for any x in the jet region; because of its location adjacent to a wall, the free jet analysis applies without any appreciable error. From the theory of conservation of momentum, the average value of the angle of expansion is 22 degrees.(2) It is (2) known also from jet studies( ) that jets discharged from a round opening form expanding cones, while jets from rectangular outlets (which is the case studied here) rapidly pass from a rectangular to an elliptical cross sectional shape at a short distance from the outlet face, and then to a circular shape, at a rate depending primarily on the jet dimension. In other words, there are three characteristic stages measured along the axis of the jet: a short preliminary stage during which the core of primary air is dissolved by mixing with room air, a second short stage during which the stream becomes circular in cross section and a third stage common to all streams whatever the shape of the air inlet is in which mixing and expansion occurs as a cone shaped stream. The first and second' stages are also short in comparison with the third stage. For engineering purposes it can be assumed that there is only one region, a semi-cone shaped stream region. Figure 2.6 shows the three phases. The velocity profile at any cross section of a free jet can be approximated by a normal probability distribution curve of the type shown in Figure 2.7.

-21where vx = average velocity of the jet at a cross section distance x from the inlet in ft./min. v = center velocity of jet at a distance x in ft./min. Xmax However, due to the existence of the wall beside the air register, the velocity profile for the stream will not be exactly half a normal probability curve. The velocity of the air at the surface of the wall will be practically zero as shown in Figure 2.8; but for engineering purposes, the normal probability curve can be used in order to calculate average temperatures over the jet cross section without an appreciable error. Experimentally it was found that:(2) -xmax. ao \ VO x and v < X max < 2.4 = - 3.2 vx - 2.8 where vo = average velocity of jet at inlet ft./min. = area of inlet in ft2 (see Figure 2.7). 2 a = discharge factor and its practical range is from 5 to 7 Define now the following: Ti(t) = average temperature of air at inlet of jet region Tx(t) = average temperature of air at a distance x Tl(t) = average temperature of air when leaving the jet region to region 2

-22T3(t) = average temperature of air in region 3 and the room initially Q = volume flow of air in ft3/min. at inlet of the jet region i = volume flow of air in ft3/min. at a distance x Q1 = volume flow of air in ft3/min. when leaving the jet region to region 2 p = density of air in lbs./ft3 Cp = specific heat of air at constant pressure in BTU/(lb) (F) Tx = time taken by a unit volume of air to travel from the inlet for a distance x in minutes T1 = time taken by a unit volume of air to travel through the whole jet region in minutes The total rate of flow of air in the stream in relation to the primary air flow issuing from the inlet can be deduced from the principle of 2 2 conservation of momentum that is: vo Ao = v Ax where the zero subscripts refer to the conditions at the nozzle and subscript x to any point measured along the axis. Since the air flow Q (volume/time) is equal to Av, one may write: x vo oEx Qo VX where E is called the entrainment ratio. It is the ratio of air flow Ex a vA -If an average value of 6 for a is used, then:

-232.8x - 6 Ao As a first step here, the average temperature at a distance x will be calculated as a function of x. If it is assumed that the air is coming in the room with a temperature Ti, the average temperature at a distance x is Tx, and that there are no losses through the inside wall then a heat balance equation can be written for the corresponding control volume to distance x. The control volume in thermodynamics is defined as a volume in space through which matter flows, and it is bounded by a control surface. The control volume here corresponds to an open system (there is matter flowing in and out of the thermodynamic system). Consider now the control volume corresponding to a distance x in the jet region shown in Figure 2.9 below. From the application of the basic first law of thermodynamics to a thermodynamic open system knowing that there is no work done in the process and neglecting the kinetic and potential energies associated with the moving air, the following heat balance equation can be written in general for any control volume in this space: Rate of change of internal energy of the control volume = Rate of change of energy associated with the mass of air flowing in - Rate of change of energy associated with the mass of air flowing out - Rate of heat transfer. For the control volume shown above in Figure 2.10, the following equation can now be written: Cp x(t)Q pC T(t )Qo +PCp( - Q)T (t) (2.1) The first term on the right hand side of Equation (2.1) represents the enthalpy flowing in from the inlet duct. The second term represents

-24Velocity Normal Probabilit Curve Actual Velocity Curve Radial Distance Figure 2.8. Velocity Distribution at a Cross-Section Normal to the Axis of the Jet. 2 Q x IXo 0 Figure 2.9. The Control Volume in the Jet Region Corresponding to a Distance x 2 Q1T1 3 T3 (Figure 2.10. C Qo)l Volume o Region 1. Figure 2.10. Control Volume of Region 1.

-25enthalpy flowing in the control volume due to the air existing in the room at temperature T. The term on the left hand side of Equation (2.1) represents the enthalpy flowing out of the control volume. There will be no heat transfer to the outside since the jet stream is facing an inside wall, and consequently there will be no appreciable rate of internal energy of this open system. Equation (2.1) can be written in the following form: QxTx(t) = QoTi(t-Tx) + (QX-Qo)T3(t) or Ti(t-Tx) 1 T (t) +~ + (1 ) T3(t) Ex Ex or Tx(t) = i -[Ti(t ) - T3(t)] + T3(t) 1.4 x Now taking the whole jet region as one control volume we have: Q1Tl(t) = QTi(t-Tl) + (Q1-Qo)T (t) (2.2) or T(t) - 4 [Ti(tTx) - T3(t)]x + T3(t) It is to be noted that the density of air and its specific heat are functions of the absolute temperature. Therefore, due to the small range of the expected variations of the absolute temperature, the air density and its specific heat are taken to be constants all over the space of the room. The average temperature of air in the jet region is seen to vary hyperbolically with the distance x as shown in Figure 2.11 below. Since the time delay in the jet region is very small compared

-26Average Air Temperature Tx Ti T1-~ —i --- - - T l 0 x x Distance along the Axis of the Jet. Figure 2.11. Variation of Tx with x.

-27with the rest of the system (which is shown in Appendix A), the equation of this region takes the following form: Q1(t) = QTi(t) + (Ql-Qo) T3(t) (2.2) The state variable diagram for the jet region follows directly as shown in Figure 2,12. 2.2.3 Region 2 Model The whole volume of region 2 will be assumed to be at a state of complete mixing, () This allows a uniform temperature assumption in this region. Call this temperature T2. From the study of air movement (2) in jets, the effect of the ceiling will not be significant until the jet approaches the ceiling to within half a jet diameter. At height H, the diameter of the jet will be 2H tan 11~ or 0.4H. Half a diameter will be 0.2H. Therefore, X will equal 0.8H Where X = length of the jet region in ft. Applying now the first law of thermodynamics to the control volume of region 2, and taking into consideration that there will be a rate of change of internal energy in this case due to the heat transfer to the outside wall we get: dT2(t) Q1PCpTl(t-T2) - Q1CpT2(t) - k2 [T2(t) - w (t) p= 2 dt (2.3) where: T2 (t) = average temperature in region 2 F~ Tw2(t) = average temperature of the inside surface of the part of the outside wall facing region 2 F~

-28T3 TiV Qo Figure 2.12. State Variable Diagram for the Jet Region. Heat Transfer to the outside 2 1 1 0.2H T2 Q1 Q T1 Figure 2.13. Control Volume of Region 2.

-29T2 time taken by a unit volume of air to travel through region 2 minutes k2 = average proportionality constant (or surface coefficient) defining the heat transfer by convection from the air in region 2 to the inside surface of the outside wall facing region 2 V2 = volume of air in region 2 (wL) (0.2H) = 0.2 wLH ft3 The first term in the left hand side of Equation (2.3) represents the energy associated with the mass of the air flowing in region 2. The second term represents the energy associated with the mass flowing out of region 2, The third term represents the energy transferred to the inside surface of the outside wall. The right hand side of the equation represents the rate of change of internal energy of the air in region 2. Again the time delay T2 of region two can be neglected with respect to the rest of the system, a result which is proved in Appendix A. Therefore: dT2(t) QPCpTl() QPCpT2(t) - k2 [T2(t) - Tw2 (t) = CpV2 d or T2(t) = alTl(t) - lT2(t) + lTw2 (t) (2.4) where Q1 Cv = __ V2 _PCpV_ C V2 CpV2 The state variable diagram for the second region follows directly from Equation (2.4).

-302.2.4 Region 3 Model Again complete mixing is expected in this region. The heat input to this region will be from region 2 by convection. The heat output will be to the inside surface of the outside wall by convection, and to the cold air duct. Applying now the first law of thermodynamics to region 3 as we did for region 2 we get: dT (t) Q1pC T (t-Tp) - QPC (t) k [T(t) - (t)] - pC V 3d 1p 5 1p3 5 5 w3 P 3 dt (2.5) where: T3 (t) = average temperature in region 3 F~ V = volume of air in region 3 ft3 3 T3 = time taken by a unit volume of air to travel through region 3 in minutes where: Tw (t) = average temperature of the inside surface of the 3 outside wall facing the third region k3 = average proportionality constant defining the heat transfer by convection to the inside surface of the outside wall facing the third region The first term in the left hand side of Equation (2.5) represents the energy associated with the mass of air flowing into region 3. The second term represents the energy associated with the mass of air flowing out of the room. The third term represents the heat transferred to the inside surface of the part of the outside wall facing region 3. The right hand side rprpesents the rate of change of internal energy of the air in region 3.

-31W21 T20 Figure 2.14. State Variable Diagram for the Second Region. Q1 / T2 Heat Transfer to the outside Q1 To Figure 2.15. Control Volume of Region 3. T1 Figure 2.16. State Variable Diagram for the Third Region.

-32By using Taylor's series expansion of T2(t-T3), and taking only the first order approximation, therefore, the state variable Equation (2.5) of the third region takes the following form: (t = 2T(t) 2T2(t)- 72T(t) -T(t) T) T (t) + 2T (t) ~2 3 ~ (2.6) where Q1 V3 Y2 - V(1 k4) Q= - Vk4 k3 k k4 = 15 QlPCp The state variable diagram for the third region follows directly in Figure 2.16. 2.3 The Wall Model In constructing the wall model, the thermal circuit techniques described in Appendix B is used. Assuming that heat flows only in one direction through the wall, the wall behaves as a distributed parameter

-33RC transmission line. As in transmission lines, it can be shown that the errors in representing the wall under the same boundary conditions varies inversely with the square of the number of T circuits used. 3) The accuracy of this model depends on the accuracy of the data of thermal conductivity and specific heat for building materials. There is some uncertainty in these data because of the porosity of the materials and uncertainty of their water content. It is apparent that the more T circuits used in describing this dynamic behavior of the wall, the higher the order of the system will be. In Chapter 5 it is shown that the order of the whole system does not affect the nature of the optimal control law but rather increases the complexity of the resulting differential equations which must be solved. It can be shown that, the number of the nonlinear differential equations to define the optimal controllers parameters is 1 + 2N + (N,)/(2((N-2).) where N is the order of the heating system. This means that a system of order three requires the solution of 10 nonlinear equations, and another of order five requires the solution of 21 nonlinear differential equations which need a large memory digital computer, With this in mind, and the fact that the state variables of the control system must be easily observed, and inexpensive to measure, it was decided to represent the wall by only one T-network. This means that a compromise between the degree of approximation and overall cost of the optimum controller is made, since this assumption puts a preliminary bound on the cost of the optimum controller. The derivation for the wall model follows: it is divided into two parts; (a) facing region 2 and (b) facing region 3.

-34Consider now Figure 2,17, which represents the thermal circuit of a one T network wall. The outside end of the wall which is exposed to the sun and wind is equivalent to a known temperature source To (the effective air temperature) and is a function of time. Ro represents the resistance whose conductance represents the heat flow between the outside surface of the wall and the surrounding atmosphere. It is a function of the surface air coefficient. The inside of the wall is facing the room heat capacity. Ri is the resistance whose conductance represents heat transfer between the air in the domestic space and the inside surface of the wall. It is also a function of the air surface coefficient. R and C are the equivalent resistance to heat transfer and heat capacity of the wall respectively. They are functions of the properties of the material of the wall. This wall system has two degrees of freedom. Usually ii, i2 (defined in the circuit of Figure 2.17) are chosen as the state variables for such circuit. However, since ii and i2 correspond to rate of heat flowing in and out of the wall and are not easily measurable in practice, another set of measurable state variables should be chosen. Let V be the temperature of the air inside the room; V1 be the temperature of the inside surface of the wall; V2 be the temperature of the outside surface of the wall. It can easily be shown that V and V1 constitute one set of state variable for this system. (7) Writing the equations for this system gives: V = ilRi + V (2.7) V1 = ilR + Vc (2.8)

V Ri V1 R Vc R V2 R Room ~4 y ~ Capacity1 2 Outsid |I" _^ I I Room Temperature Figure 2.17. The Wall Thermal Circuit.

-36c = i2R +V2 (2.9) V2 = i2Ro + Eo (2.10) t 1 (i2')dt = Vc (2.11) i-J!miLnating iI, i2, and Vc between Equations(2.7), (2.8), (2.9),;id. (2.10) yields: R R. + R R R [V1 - V+ EO ] (2.12) V2 =R +.-R R i R 0 i Again from Equations (2.7), (2.8), (2.9) and (2.11), i, i2, and Vc are eliminated to yield: B V - V- V2 _-E V = R (V - V1) + i [ - ] dt (2.13) 1 R^i C o Ri Ro Differentiating Equation (2.13) with respect to time gives: R + R. dV V- V1 V2- E R d1 +( i 1 [ (2.1) R. dt R. dt C R. R By substituting V expressed by Equation (2.12) into Equation (2.14), transposing and rearranging terms in the resulting equation gives: dV - CR(R+Ro) V + C(Ri+R)(R+Ro) - (2R+Ro)V + (2R+RoRi)Vi - RiE = 0 at dt (2.15) Using Equation (2.15), the state variable equations for parts of the wall facing regions two and three are stated in the following, where subscripts 2 and 3 refer to regions two and three respectively.

-37For region 2: In this case T2 corresponds to V and Tw2 corresponds to V1 o Therefore: ~dT2 ddT d2 w C2R2(R2 + Ro2) - + C2 (Ri + R2) - (2R2 + Ro2) T2 dt 2 + (2R2 + Ro + Ri2) T - Ri2T = (2.16) From this differential equation, the state variable diagram for the part of the wall facing region 2 follows as shown in the next section. For region 3: Again T3 corresponds to V and Tw3 corresponds to V1: dT3 w3 - cR3(R3 +Ro3) - + C3(R + R3) dt (2R3 + Ro0) T3 + (2R + Ro3 + - Ri3 Tw 0 (2.17) where: C2 = equivalent heat capacity of the part of the wall facing region 2 =?,2 equivalent resistance of the part of the wall facing region 2 -2 = equivalent outside air to surface resistance of the part of the wall facing region 2 Ri2 = equivalent inside air to surface resistance of the part of the wall facing region 2 T = outside equivalent temperature o C3 = equivalent heat capacity of the part of the wall facing region 3 R3 = equivalent resistance of the part of the wall facing region 3

— 38R3 = equivalent outside air to surface resistance of the part of the wall facing region 3 R = equivalent inside air to surface resistance of the part -i of the wall facing region 3 The state variable diagram for the part of the wall facing region 3 follows and is shown in the next section. In Sections (2.1.3), (2.1.4), and (2.1.5) it was shown that: Tl(t) = l Ti(t) + (1 - ) T (t) (2.2) 2(t) = al Tl(t) - 1T2(t) + ZTw (t) (2.4) T(t) = - 2Tl(t) + 2T2(t) - 2T3(t) - 2T (t) (26) + t T (t) 2wv By substituting T2 in Equation (2.4) into Equation (2.16), transposing and rearranging gives: Tw2 = -Tw2 T1 T- 2 + 30To (2.18) Similarly for iw by using Equations (2.6) and (2.17), therefore: T T3 + 4T1 2 843 = 4T2 + 4to (2 19) 4w Y4T2 4T3 4T T w3 - % 3 2 where: 3 c= ^(RI2+R2) [71C2I)(R 2 + Ro2) + 2R2 + Ro2 + Ri2]

-39alC2R2(R2(2 + Ro2) C2(Ri2 + R2) 3 c2(R2- ) [c22(R2+ Ro2) + (2R2+ Ro2)] R i2 C2(Ri2 + R2) C64 = )[2C3R3(R3+ Ro) + (2R3+Ro3+Ri3)] c423 13(3 o3) C3(Ri3 + R3) n^~1 _, -4~ ~ [2c3R3(R3+ R+)] CQ(R +R5) 2.4 The Room System Model Summarizing the set of equations representing the room space and walls is: Q Q T, = -2R + (1 - - ) T (2.2) i2 = a - 1T2 + 71T (2.4)

-40T2 Tw Figure 2.18. State Variable Diagram for Part of the Wall Facing Region 2. 4 l T 2, _____,^w To Figure 2.19. State Variable Diagram for Part of the Wall Facing Region 3.

-41T3 = 2T1 + 2T2 - 72T3 - 62Tw + 52Tw (2.6) Tw = - T T1 - 3T2T + To (2.18) w3 = -Tw + 4T1 Y4T2 -4 - 4Tw2+ 4T (2.19) This is a set of linear first order differential equations with constant coefficients. This set represents the room system, Ti is the input and T1, T or T may be considered as the output of the system. Using these differential equations the state variable diagram for the room system follows as shown in Figure 2.20. This mathematical representation of the room is the outstanding problem in the theoretical treatment of the operation of a warmed air domestic space heating. The adequacy of the model will be checked experimentally in Chapter III. The next chapters will also describe the incorporation of this model in an analog computer model of the complete heating system, to investigate the dependance of the performance of the heating system on the system parameters.

^] I t ~ I. QiTi~~~~~~~~~~~~~~~~~~~~~~~ 2o~ ~~~~~~~ro Figure 2.20. State Variable Diagram for the Room System.

CHAPTER III THE EXPERIMENTAL CHECK OF THE ROOM MODEL This chapter is devoted to the description of the complete experimental work done towards checking the validity of the mathematical model of the domestic space. A thorough evaluation of the applicability of the mathematical model of the room requires a parallel comparison between a physical installation and the mathematical model, To determine the validity of the mathematical model predictions, measurements of the thermal environment, and the response of a test room were made during February, 1966. The analoge computer results compare favorably with the observed data. Graphical comparisons of results are included. 3.1 System Investigated 301.1 General Description The experimental room model consisted of a test room located in the warehouse portion of the Maxitrol Company, Southfield, Michigan. Instruments and controls for obtaining data, electric heaters and air conditioner for the room air was an integral part of the installations. A view of the interior of the test room looking south, is shown in Figure 3.1. The room was 12 ft. wide and 12 ft. deep and had a 10 ft. height. These were all outside dimensions. The north wall, not shown, was a duplicate of the south wall with the exception of a door of 3 ft. width and 6.5 ft. height centered in that wall. Of course, neglecting any leakage of air through this door, the south and north walls can be considered identical since the door construction is the same as the wall construction, all walls are of similar details, and as seen in Figure 3.1 -43

-44A closed Air Inlet ~ An open Air Inlet \^\ P J~IControl Observation Windows _go~ut Wall Air Returns through Walls Figure 3.1. View of Interior of the Maxitrol Room. 12' Observation, windows \ I ^ ~Control Wall Window /I L_.1 I / 9,67 - Figure 3,2. The Maxitrol Test Roori with Bimens ions shown

-45the south wall contains neither windows nor doors. The west wall (to the right in Figure 3.1) had only a control window 1' x 1' which connects between inside sensing equipment and the outside recording equipment. This window is also of the same construction as the walls, with the exception of holes through which cables can be run. This window opens and closes like a door. Of course, it was kept closed all the time during the test. The east wall (to the left in Figure 3 1) had two small glass windows each 1/2' x 1/2' to permit observation of the interior from the outside. This room was designed for the Maxitrol purposes to have a special ceiling with holes in it. For the purpose of this test, this ceiling was removed, and the east air duct entrance was closed by a plate of metal so that only one jet form of air would be studied. This room was designed also to have the air return through the walls as shown in Figure 3.1. In otherwords the air constituted a layer in the construction design of the walls. The air which was introduced into the room was conditioned by equipment which consisted of electric heating coils for heating the air and a two ton refrigerating unit, utilizing chilled water coils, for cooling the air. A variable speed fan capable of delivering flow rates of about 3000 ft3/m., returned the room air from a return intake through ductwork to the conditioner. After the air passed through the conditioner the fan returned the air to the room through supply ductwork and the outlet arrangement. The electric heaters and the air conditioner system was located outside beside the south wall. The control and recording equipments were placed outside the west wall.

-463.1.2 Construction Details All the walls and the door were made of: 3/8" exterior plywood, 1" fiber glass, 0.025" sheet of aluminum, 1-1/4" air space, 0,032" sheet of aluminum, and 1/8" tempered hardboard (tempered masonite). The ceiling consisted of: 3/8" exterior plywood, 1" fiberglass, 0.032" aluminum sheet, 3" rolled fiber glass insulation, 0.032" aluminum sheet, and 1/8" tempered hardboard. The floor was made of: 1-1/2" exterior plywood, 4" air, and 3-5/8" wood. These construction details are summarized in Figure 3.3. 3.1.3 Instrumentation The test room was provided with instrumentation for measuring the physical variables of the system that have an important influence on the thermal behavior. Measurements were made of the following variables: 1. Inside and outside temperatures 2. Inside and outside air velocities Measurements were recorded continuously during the test period. 3.1.3.1 Temperature Measurements Thermistors of the type RL25J1T were used to measure the temperature. The thermistors were used in series with a 3 x 105 ohms linear resistance in order to avoid selfheating. The electronic recorders were calibrated to read the voltage drop across the thermistor, and from the calibration chart of the thermistor the corresponding values of temperatures were determined. Figure 3.4 shows the setup used to measure the temperature.

-47h | |\ 3/8" Plywod e _ ~A\V I I^ 1" Fiberglass to _ ~ \\ 0.025" Sheet of Aluminum X _- ~~~- - \ 1-1/4" Air Space 00,032" Sheet of Aluminum ^_ _. ~~~ _^ \1/8" Tempered Hardboard (Tempered Masonite) a. Wall Construction \\\\3/8" Plywood - ~\\\l1" Fiberglass Inside Air \\\ 0.032" Aluminum Sheet " Rolled Fiber Glass Insulation 0.032" Aluminum Sheet \/8" Tempered Hardboard b. Ceiling Construction Inside Air~ 3-5/8" Wood'"1 4?" Air 1.5" Plywood c. Floor Construction Figure 3.3. Construction Details, of the Test Room Boundaries,

-483 x10 5 E LECTRONIC RECORDER THERMISTER _ 10 v Figure 3.4. The Setup for Temperature Measurement.

-493.1.3.2 Velocity and Rate of Flow Measurements A velometer of the type Alnor was used to measure the velocity and rate of flow of air at different points in the room. 3.2 Calculation of System Parameters 3.2.1 Specific Heat of Air (Cp) The specific heat of air in the region from 20~F to 120"F remains constant at 0.24 Btu/lb~F. Therefore: Cp = 0.24 Btu/lbF. 3.2.2 Density of Air (p) It varies slightly with temperature in the range between 20~F and 120"F. An average value over this range is taken and is equal to 0.0749 lb./ft. 3.2.3 Volume Calculations According to the mathematical model described in Chapter II the test room is drawn schematically in Figure 3.5 below, showing the corresponding regions and their dimensions. Using these dimensions the volumes of the different regions are calculated. Details are in Appendix D. V1 = Volume of region 1 = 30.7 ft3 V2 = Volume of region 2 = 270 ft3 V3 = Volume of region 3 = 1049.3 ft3

-5022 93 93 233" 23.3 a, Elevation View of the Test Room. 11.76' 11.76' b. Plan View of the Test Room. Figure 3.5. The Test Room with Regions Shown (Dimension indicated are all internal).

-513.2.4 Thermal Conductivities Details of calculations are in Appendix D, where it is found that: kw = equivalent thermal conductivity of the wall = 0.0238 Btu/(ft. OF. hr.) kc = equivalent thermal conductivity of the ceiling 0.0267 Btu/(ft. OF. hr.) kf = equivalent thermal conductivity of the floor = 0.0338 Btu/(ft. F. hr.) 3.2.5 Thermal Resistances The following parameters were calculated in Appendix D: R = 1.975 OF/(Btu/min.), 2 Ro2 = 17.8 ~F/(Btu/min.), Ri = 17.8 ~F/(Btu/min.), R = 0.66 ~F/(Btu/min.), RH = 8.2 OF/(Btu/min.), Ri 8.2 ~F/(Btu/min.). i2 3.2.6 Equivalent Capacity of the Room Boundaries As explained in Chapter II an equivalent capacity for the walls and the ceiling and the floor has to be estinated by some means. Like the electrical capacity, the heat capacity for a small piece of material is defined as the ratio of the quantity of heat stored in this material over an interval of time to its temperature rise during that period of time. This definition is used here to obtain experimentally an estimate for the capacitance of the room boundaries.

-52The quantity of heat added over an interval of time = heat added to the room through inlet duct - heat leaving the room outlet duct - heat stored in the room due to its capacity - heat lost to the outside during that time. The average temperature rise of the mass of the walls, ceiling and floor will be estimated by taking the average rise of temperature of the inside and outside surface temperatures. A test was made on the room model to determine an estimate for the equivalent capacity of the building. It consisted of the following: Heated air was pumped in at the rate of 125 ft3/min. Measurements were taken for: the input temperature to the room Ti, the average inside surface temperature and the average outside surface temperature and the average outside surface temperature of the walls. Plots of these measurements are shown in Figures 3.6, 3.7, 3.8, and 3.9. The outside air temperature was held constant at 73~F. To calculate the heat input to the system, consider the interval of time starting from the time 20 minutes to 40 minutes. The heat input starting from the time 20 minutes to 40 minutes corresponds to the area under the input temperature curve shown in Figure 3 6. Heat input during 20 minutes = pCpQo x Area under the Ti curve = 0.0749 x 024 x 125 x 1600 = 3590 Btu The heat carried out in the return duct which is placed in region 3 corresponds to the area under the curve of Figure.7 between the time 20 minutes and 40 minutes.

-5390 ~ ~~C~ ~1~y 90 80 80 50-. T50 TIME i MIN. TIME, MIN. Figure 3.6. Input Temperature for Figure 3.7. Average Temperature of Capacity Measurement, Region 3 for Capacity Measurement. Iw w w, w w 70 ~~ 70 w w La.. 60 60 w w 50 5 0 10 20 30 40 50 60 0 10 20 30 40 50 60 TIME, MIN. TIME, MIN. Figure 3.6. Insidput Temper aturface forWall Te- Figure 3.7. OutsAverage Temperature ofeperature for Capacity Measurement. Region 3 for Capacity ~~~Measurement,~~ Measurement. Z v, Measurement. Measurement.

Heat leaving the room during 20 minutes = pC Q x Area under the T1 curve = 0.0749 x 0.24 x 125 x 1400 = 3142 Btu Heat stored in region 3 = pCpV AT P 3 3 = 0.0749 x 0.24 x 1049.3 x 6.5 = 123 Btu where AT = increase of average temperature T3 of the third region over the 20 minutes interval. Heat stored in region 2 = pCpV2 AT2 = 0.0749 x 0.24 x 270 x 6 = 29 Btu where AT2 = average increase of temperature T2 of the second region over the 20 minutes interval. Similarly heat stored in region 1 = 4.33 Btu Heat transfered to the outside k2 (Toutside - T) x time surface + k3 (Toutside - To) x time surface = 3.93 Btu where To = outside air temperature = 73~F Heat stored in the room boundaries (capacitor) = 3590 - 3142 - 123 - 29 - 4.33 - 3.93 = 288 Btu

-55Average temperature rise of the equivalent capacitor = 8~F As defined before the room boundaries equivalent capacity is: _288 Btu = 36.2 8 ~F. This equivalent capacity is divided in proportion to the areas facing region 2 and region 3 to give C2 and C3 respectively. Therefore: C2 = 225 x 36.2 = 115 2 723 0F. 3 = 455x 36.2 = 22.8 Btu 723 ~F. This estimate of the capacity is checked using another experimental set-up in Appendix C. The system parameters are calculated using the proporties of the individual layers. An electrical analog is built and the equivalent capacity is estimated. 3.3 The Experimental Test The main purpose of the test is to check the validity of the mathematical model of the room developed in Chapter II. The room model developed was shown to be linear and time invariant in the range of domestic air temperatures. This followed from the fact that the air parameters are practically constant over this range of operating temperatures. They are independent of time or temperature throughout the range of domestic air temperatures. It also followed from the lumped parameter approximation of the room boundaries, and the assumption that these parameters are independent of time and temperature. To check

-56the validity of this model over the operating range of temperatures an input of temperature in this operating range is applied to the system and the output temperature response is observed. Since it is desirable to use a level of temperature that is used in actual heating systems a 40-110~F temperature range input was used in the test. A comparison then is made between the output of the system observed experimentally and the output from the corresponding analytical model. In Section 3.4 the results of a computer simulation is given for the corresponding analytical model. In this section the description of the test and its results are presented. The input temperature Ti was obtained by turning the heater on and the refrigerator off at the same time. The temperatures at the input, in region 1, in region 2, in region 3 and out side the domestic space were connected to the recording equipments. Figure 3.10 shows the location of the points where temperatures were measured. The state of the system was initially at: T = 53.80F 20 T30 = 54.1~F T = 57~F w30 and was in steady state. This state was created by having the refrigerating unit on, the heaters off, and ajusting the blower to give a constant rate of air flow of about 1340 ft3/min. Starting with these initial conditions, the heaters were turned on and the refrigerator off. The air flow was kept constant at 1340 ft3/min. Recordings for temperatures in

47-b / IA\,^/ \ i /A1 \ /.B \ / GA I eE ~ eE,.. /DOF.F D Side view Elevation view *F Ie Be [)i,A Ge'C eE eH eD Plan view Figure 3.10. Location of Points at which Temperatures were Taken During the Test.

-58the different regions were continuously taken. Three points in every region were observed: TA, TB, and TC in region 1: TD, TE, and TF in region 2; TG TH, and TI in region 3. The input temperature Ti was also taken. The time delays in the different regions were recorded. The outside air temperature To was essentially kept constant at 73~Fo The outside wind speed was basically zero, since the machine shop was closed at all times. The test ran for about 120 minutes, at which steady state conditions were established. Plots of these recordings are shown in Appendix Do 3.4 Computer Simulation The purpose of this section is to simulate the same conditions of the experiment on the analog computer set up for the mathematical model of the domestic room developed in Chapter II, and compare the experimental results versus the analytical results in the transient and steady state conditions, Before the electronic analog computer simulation, it is necessary first to calculate the parameters of the Maxitrol room model i.e. the parameters in Equations (2.2), (2.4), (2.6), (2.18), and (2.19) which corresponds to the particular experimental model. The calculations of these parameters are straightforward, and they follow directly from the results of Section 3.2. 01 = 2. 25 2 = 0.769 03 = 0.0112 04 = 0.0107 P1 = 2.26 Y2 = 0.584 83 = 0.225 4 = 0.0142 Y = 0.0116 62 = 0.00098 y3 = 0.231 74 = 0.0518 a~ = 0. 191 2 = oo00646 6 = 0.00452 = 0.0435

-5954 = o =.39 14 = 0.0047 1 - = 0.61 Therefore Equations (2.2), (2.4), (2.6), (2.18), and (2.19) take the following form for this experimental model: T = 0.39 T3 + 0.61 Ti T2 = 225 T1 - 2.26 T2 + 0.0116 T2 T = 0.769 T -0.191T1 - 0.584 T- 000098 Tw + 0.00646 Tw t2 = - 0.225 T1 - 0.231 T2 - 0.0112 Tw - 0.00452 To Tw = - 00.0518 T + 0.0142 T 0.0435 T - 0.0107 Tw + 0.0047 To w' 2 1 3... The unit used for time in the equation above is the minute, i.e. if these equations were to be simulated on the analog computer then every second on the computer represents one minute of real time. However, a suitable scale for this problem is 1:240 i.e. 1 second of the computer time corresponds to 4 minutes of the real time. The scale for the temperature is 0.5 volts corresponds 1~F. This is accomplished in this setup by scaling only all inputs and initial conditions. The system equations can be written after scaling in the following form: T1 = 0.39 T3 + 0.61 Ti (3.1) jT = 9 T - 9.04 T2 + 0.0464 T (3.2) T~~~~~~~~~~~~~32)1

-60N M 0 K cJ'K L Figure 3.11. Location of Points on the Inside Surface of the South Wall, Where Temperatures were Recorded During the Test.

-61T3 = 3.076 T2 - 0.764 T - 2.336 T - 0.00392 T2 + 0.02584 Tw3 5 2 1 5 w2 w5 (3.3) Tw2 = - 0.9 T2 - 0,926 T - 0.04504 Tw2+ 0.01808 To (3.5) T = - 0.2072 2 + 0.0568 T1 - 0.174T- 0.043 Tw + 0.0189 To (3.6) To simulate the same experimental conditions on the analog computer, the generation of the input function Ti similar to that applied in the experiment is necessary in order to compare the outputs. The actual input obtained experimentally is shown in Figure 3.13. Using this plot a piecewise linearization is obtained and shown in Figure 3.12. A computer diode function generator is then used to generate such a piecewise linearization of the input function T.. The computer set-up for this purpose is shown in Figure 3.14, and the function generated by this setup is shown in Figure 3.13. T;i'F~ 105 100 831 — - 66.5 46.25I Real time!i i$___________!____I I_ 1____in m inutes 02 5 13 51 93 Figure 3.12. Piecewise Linearization of the Input Function.

110 100 90 -- ~ ~~ ~~~~~ 90 The solid line represents Ti generated -L 80 ~ i on the analog computer. The points re_' ~~~~~f ~present the experimental values. 50 ~ I I I I I I I II I I I I U 70 60 50 40 0 10 20 30 40 50 60 70 80 90 100 110 120 TIME, MIN. Figure 3.13. The Input Temperature.

DFG -100v T; -l G4.25 \\ \ 1 GK F6 Fl'~ ^+X ^ / / ^.231) F2 -<-+00v Figure 3.14. Analog Computer Diagram to Generate the Function T. DFG: Diode Function Generator.

-64Using now the Equations (3.1), (3.2), (3.3), (3.4), and (3.5) the analog computer set-up for the experimental room model is constructed and shown in Figure 3.15. The description of the computer set-up here is unnecessary, since it follows standard practice. TABLE 3. COMPUTER POTENTIOMETER SETTINGS OF FIGURE 3.16 A6 0.9 B5 0.234 C6 0.265 1 Gl 0.27 A7 0.046 B6 0.026 D1 0.017 G2 0.926 B 0. 904 B7 0. 004 D2 0.043 G3 0.900 B2 0.269 C1 0.057 F1 0.25 G4 0.007 B3 0.308 C2 0.207 I F2 0.231 G5 0.045 B4 0.764 C3 0.007 F6 0.04 i G6 0.265 Computer runs were made on the Applied Dynamics analog computer for the generated input function Ti. The outputs of the system T1, T2 and T3 were recorded using electronic recorders. Plots of these recordings are shown in Figures 3.17, 3.18, and 3.19 by the dotted curves. 3..5 Evaluation of Model Examination of Figure. 3.16 indicates that the average temperature response T1 in region 1 experimentally and analytically shows the steady state error to be 4.02% and the maximum error during the transient period to be 3.53%. Figure 3.17 which also shows the average temperature in the second region, indicates a steady state error of 4.07% and a maximum transient error of 7.46%. Similarly, from Figure 3.18, the average

A5 ~0~~~~~~~~~~~~~~~6 B4 Bl B5:i4 42~ A4 A AG ~T126 T2 -T 2 1 Ti +Tl a ~~~~-T3+ 1B 021C.6 2 B 2\.2 Gl A7 +100v ~ (.017) Dl G2 0 + B6 Figure 3.1$. Analog Computer Diagram for the Room System. C 1.057 G3 +10OV, 2 2 c6 (25)c +100v 1 100V 0 5D Figure 3.15. Analog Computer Diagram for the Room System.

100 -o- -0 900 8 0 90 80 20_30 40 50 60_70_80_90__00____ __20_130 80 ~~70t7: Figure 6 AergTmet -— 1 Experimental Analytical 60 - 5010 20 30 40 50 60 70 80 90 100 110 120 130 TIME, MIN. Figure 3.16. Average Temperature in Region 1

-67100 ~~~ 501 I00 II I I L 00 0 90 ~ " 0" / Figure 3,17. Avr —- Experimental I -- Analytical 60....... ~ ~ ~ ~ ~ ~ ~ ~ 50 ~ 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 TIME, MIN. Figure 3.17. Average Temperature in Region 2. 90 _0 1 70 ~ ~^^~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 90 -- - Experimental - Analytical 60 0 10 20 30 40 50 60 70 90 100 110 120130140 TIME, MIN. Figure 3.18. Average Temperature in Region 3.

-68Li. 80 50 ~ ~ - _ ~ ~ 40- 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 TIME, MIN. Figure 3.19. Average Wall Surface Temperature Facing Region 2. 80 --- Experimental - Analytical 40 ~~_ 10 20 30 40 50 60 70 80 90 100 i10 120 130 140 150 TIME, MIN. Figure 3.20. Average Wall Surface Temperature Facing Region 3.

-69temperature in the third region has a steady state error of 1.48% and a maximum transient error of 3.1%. Figures 3.19 and 3.20 display the corresponding graphs for the wall surface temperatures T, T w2 w3 The steady state errors are 3.9% and 4.66% respectively. The maximum error during the transient period is 5% for Tw and 4.23% for T. The analytical temperature response when compared to the observed temperature response of the test proves, to within a reasonable average error, the validity of the mathematical model, for each individual temperature response. Due to unavoidable errors, such as instrumention errors, and errors resulting from the difficulty in obtaining the exact characteristics of the building material, there is a deviation of results from the physical situation to the analytical model. It is estimated (work done by Buchberg at UCLA)(9'10) that all these errors pITt together can easily result in errors of the order of five percent in any temperature quantity computed. These results, however, in most temperature responses, indicate a percentage error lower than 5%, with the exception of one transient case, where the error reached 7,46% at one instant. The steady state responses were well below 5%. The results also make apparent that for each analytical case as well as for the experimental case, T1 is greater than T2, T2 is greater than T3, and T3 is higher than both T2 and T. However, the differences between T1, T2, and T3 are comparable within 5%. In addition, the difference between T and T is also comparable within 5%. W2 W3 One objective of this dissertation, however, is to study the room temperature variations of the conventional heating system, and compare it with that of the optimal heating system. The jet region, due to

-70its small volume is of no practical importance here. As a matter of fact, region 3 occupies the greater portion of the room area, and is the region of practical interest. For this type of study, the mathematical model seems to be satisfactory, since in these situations differentiation between T2 and T3 is not of major importance. However, it appears that the overall mathematical model of the heating system is greatly simplified if the temperatures T1, T2, and T3 are considered to be one uniform temperature, TR, and if the temperatures Tw and Tw are considered w23 to be one temperature, Tw. This approximation, a result of the experimental data, is of great engineering significance, as it reduces the number of nonlinear equations representing the parameters of the optimal law from twenty-one to ten. It also reduces the number of state vraiables to be measured. In other words, such approximation also represents a preliminary bound on the cost of the optimum heating system. Having this in mind, as well as results of the experimental data, it is decided to combine T1, T2, and T3, and were treated as one uniform temperature, TR, and treat T, and T3 as one uniform temperature, Tw The approximate mathematical model of the room now assumes the following form: pCpV r1R = (pCpQ Ti - pCpQ TR) - k(TR - Tw) (3.1) - CR(R + Ro) R + C(R + R) T - (2R + R) TR + (2R + R + Ri) T - Ri T 0 (3.2) where V = volume of the room Q = rate of flow of air

-71k = average proportionality constant defining the heat transfer by convection to the inside surface of the outside wall C = equivalent heat capacity of the outside wall R = equivalent resistance of the outside wall Ro = equivalent outside air to surface resistance of the outside wall Ri = equivalent inside air to surface resistance of the outside wall Equation (3.1) represents the simplified heat balance of the room space when T = 2= T = TR. Equation (3.2) represents the mathematical model of the wall when Tw = Tw = Tw as derived in Section 2.3.

CHAPTER IV THE CONVENTIONAL HEATING SYSTEM ANALYSIS In this chapter the analytical model of the room together with the mathematical model of the components of a conventional heating system is used to analyze the conventional heating system. The purpose of this analysis is to study the performance of the conventional heating system using the analytical model, and also to be able to compare it with any optimal or sub-optimal system to be devised later on. Mathematical models for the components of the heating system are developed. The models of the elements thus derived together with the room mathematical model are then to be simulated on an electronic differential analyzer. Typical heating system parameter variations are studied and analyzed using the validated domestic room model. 4.1 The Mathematical Model of the Controller The controller's function is to maintain the enclosure air temperature at a desired level by controlling the heat energy input to the heating system from the heat exchanger in an on-off fashion, The class of controllers used in most present day domestic heating systems consists of thermostats which are of the bimetallic element type with and without compensation. These thermostats yield a mechanical output variable which is used as the primary feedback element in many warm air heating systems currently in existence, Since the rate of change of temperature of the bimetallic element is proportional to the heat flowing into it (follows from the first law of thermodynamics), it can be stated that the rate of change -72

-73of heat of the bimetallic element is proportional to the difference in temperature between the thermostat element temperature and the temperature of the surroundings. Therefore: mc dT(t) UA [Ts(t) - Tm(t)] (4.1) dt where m = mass of thermostat element c = specific heat of thermostat element A = boundary area normal to heat flow U = overall heat transfer coefficient Tm = thermostat element temperature Ts = average temperature surrounding the thermostat element. Since a linear relationship exists between the mechanical movement of the thermostat element and its temperature Tm(t), then: xl(t) k1 Tm(t) (4.2) where kl = proportionality constant xl(t) = represents the mechanical displacement of the thermostat element from some suitable position Similarly if a linear relation exists between the set point temperature and contact position then the contact position as determined by the set point temperature may also be expressed in a form similar to Equation (4.2), above. In other words x2(t) = k2 TI(t) (4.3)

-74where x2(t) = contact position k2 = constant of proportionality TI(t) = set point temperature. This relationship may be visualized by assuming that the set point positions the contact by means of a screw-type adjustment as illustrated in Figure 4l1. Combining Equations (4.2), and (4.3) the mechanical output x(t) = x2(t) - xl(t) is given by x(t) = k2 T(t) - kl Tm(t) (4.4) Using Equations (4.1), (4.2), (4.3), and (4.4), the state variable diagram for the thermostat is derived and is shown in Figure 4.2. TT in the diagram is the thermostat time constant and is defined by: mc 7T = iUA Included in this diagram is the correction signal produced by the on-off contact or which is +1 if the furnace is on and 0 if it is off. This is shown in greater detail in Figure 4.3 where h represents the hystersis range. This means the furnace will be on (+1) until the room temperature exceeds a certain value, then the furnace turns off (0) until temperature drops to a certain value at which time the cycle repeats itself. 4,2 The Mathematical Model of the Heat Exchanger There are different kinds of gas furnaces in practice, A typical gas furnace that is most commonly used would be the vertical tube

-75Voltage Supply, V, Gas Valve x 1 //////\7 )/ 7 —7 77 Reference Figure 4.1. A Bimetallic Thermostat. T TI Figure 4.2. State Variable Diagram of the Thermostat. e(t) +1 h Figure 4.3. Hysteresis Function.

-76combustion chamber shown in Figure 4.4. Figure 4.5 shows a heat transfer schematic for this type of furnace, Considering now the thermodynamic control volumes to be the material of the heat exchanger as one and the air circulating around the exchanger as another, the following set of equations can be written for the system by the help of the first law of thermodynamics: If(t) = hf D Y(Tf(t) - T (t)) (4.5) Ie(t) = ha De Y(Te(t) - Ta(t)) (4.6) Ic(t) = PCpQoTR(t) (4.7) a(t) = PaCaQoTa(t) (4.8) d Te(t) If(t) - Ie(t) + wce (4.9) dt d Ta(t) 1 ~d T= - [Ie(t) - Id(t) + Ic(t)] (4.10) dt wa a where: De = Diameter of the heat exchanger Y = length of the heat exchanger Tf(t) = average temperature of flame Te(t) = average temperature of the heat exchanger wall Ta(t) = average temperature of the cold air If(t) = heattransferred from flue gas to exchanger wall Ie(t) = heat transferred from exchanger wall to air Ic(t) = heat coming in with the cold air Ia(t) = heat carried away to the hot air duct pa = average air density in the heat exchanger hf = heat exchanger coefficient between flame and heat exchanger material ha = heat transfer coefficient between heat exchanger material and circulating air.

-77 - Insulat ion~ Vertical Metallic Tube Burner Gas Air Mixturer from Air from secondary air - l cold air duct Figure 4.4. A Simple Diagram for a Typical Gas Furnace.

-78Ia Ta Te Tf ITf $ Ie If aS Ic Figure 4.5. Heat Transfer in the Furnace.

-79Ca = average specific heat of air in the exchanger QO = average flow rate of air in furnace. This set of Equation (4.5) to (4.10) may be summarized by the following two differential equations: d Te(t) hf $DeY[Tf(t)-Te(t)] = ha tDeY[Te(t)-Ta(t)] + WeCe (4.11) dt a dt ha DeY[Te(t)-Ta(t)] PaCaQoTa(t) + PaCpQoTR(t) (4e12) Transposing and rearranging these equations may be expressed in the following form: Te(t) = a1 Te(t) + a T(t) + a3 Tf(t) (4.13) Ta(t) = a4 Te(t) + a5 Ta(t) + a6 TR(t) (4.14) For more details regarding the furnace model the reader is referred to work done by: L.F. Kazda, Graham Casserly, Adel Eltimsahy, and Ronald Spooner (12) In practice, however, it was found that the heat capacity of the air flowing around the exchanger wall is very small compared with the capacity of the heat exchanger. (12) It is a fact that the time constant of the heat exchanger (which is defined as the time necessary for the air in the heat exchanger to rise to 63.2% of its final value when the flame temperature changes abruptly)is proportional to its capacity, Therefore the heat exchanger material capacity contributes to most of the exchanger time constant. It is then justified from the engineering point of view to use the heat exchanger material capacity

TR a1 Te o Figure 4.6. State Variable Diagram for the Heat Exchanger.

.81only for the purpose of studying the effect of varying the time constant of the furnace. Therefore the state variable equations for the simplified model take the form: T(t) = a1 T(t) + a2 Ta(t) + a3 Tf(t) (4.15) 0 = a7 Te(t) + Ta(t) + a8TR(t) (4.16) where: a _ w_ [hf -jDeY + ha TtDeY] haiDeY a2 -- T a^a e 2 WeCe hfrDeY w c e e h, JTDeY a = _ ___a e____ 7 PacaQo + ha iDeY a8 =- C- ~ PaCaQo + ha iDeY From these equations the heat exchanger state variable diagram shown in Figure 4.6 below is obtained. 4.3 The Mathematical Model of the Hot Air Duct The warm air heating and air conditioning systems ducts are the vehicle by which air is directed from a primary source to the various parts of the house and, where desired, returned to the primary source, In accomplishing this job the ducts need only to contain and direct the air with minimum heat loss.

-82TR Ta Ti I I Figure 4.7. Schematic of an Air Duct. If the temperature in the remainder of the house is assumed to be at TR, the following imperical formula(36) (ASHRAE 1963) which represents the duct mathematically can be written as: ) Ti(t)(r+l) - 2TR(t) Ta(t) ~ (r-1) where 28.8 Ad Vp r = - A for rectangular ducts UpI 7.2 doVp r =. 0 P for round ducts UO Ad = cross sectional area of duct ft2 a p = density of air lbs./ft3 do = Diameter of round duct, ft U = overall coefficient of heat transfer for the duct wall Btu per (hr) (-ft2)(F~) & = length of duct, ft Ta(t) = average temperature of air when leaving the heat exchanger. p = perimeter of the duct

-83 Therefore the state variable diagram for the hot air duct follows as shown in Figure 4.8. 4.4 The Mathematical Model of the Gas Control Valve The control valve which releases the gas to the heat exchanger will in general be operated by a voltage6 el(t). It will possess a coil containing resistance R and inductance L. The valve will contain mass M, viscous damping a and a restoring spring K. The valve gas line pressure p is exposed to a surface area Ag. Applying Lagrange's equations to the system yields a set of nonlinear differential equations. When linearized about a set point operating value, the following differential equations are obtained: ei(t) = Ri + L Kv dt dt p(t)Ag = M 2 + a + (K-Kv) x + Kvi dt2 dt It is true however, that the dynamics of the gas control valve can be neglected when placed in the conventional heating system. 4,5 Analog Computer Model Having the system all formulated in the state variable form the next step is to construct a suitable analog computer setup in order to analyze the heating system. Provision is made in order to study the dependence of the performance of the heating system on the system main parameters: thermostat time constant, thermostat hysteresis, and heat exchanger time constant. The procedure of building an analog computer for each component follows standard rules. Figure 4.9 represents an

^~\^ ~ ^ i Tigare h-.Q. S-fcate Variable Diagram for -the Ho-b Air Duct. T3 Figure 4.8. State Variable Diagram f or the Hot Air Duct.

-TR K.1 ( -T D3 El D4 E2 Figure 4.9. Analog Computer Diagram for the Thermostat. TR TR E6 1 1 1J1 -Tf ~]\T~ \ ^Tj I 2 ~ Fl ~(.98 F'~ c4 E7 xF4 E5 Figure 4.10. Analog Computer Diagram for the Heat Exchanger.

-86analog computer diagram for an on-off thermostat. The hysteresis is varied by potentiometer E2. The time constant is controlled by using potentiometer D3 and D4. Potentiometer El controlls the gain of the thermostat. For a typical heat exchanger: hf = 0.0667 Btu/(min) (ft2) (OF) ha = 0.111 Btu/(min) (ft2) (~F) Ca = 0.241 Btu/(lb) (F) pa = 0.064 lbs/ft3 density of heat = 475 lbs/ft exc hqn.er material..-X =- 0.11. Btu/lb/OF Using these parameters, an analog computer diagram for the heat exchanger is constructed as shown in Figure 4.10. Provision is made to vary the time constant of the furnace. For the hot air duct, choose the following: Ad = 1ft2, p = 4 ft, p = 0.064 lbs/ft3 & = 20 ft, v = 1340 ft/min., U = 0.49 Btu/(hr) (ft) (~F) r = 63 The computer set up for the hot air duct follows as shown in Figure 4.10 with the heat exchanger. The analog computer set up for the simplified room is shown in Figure 4.11. 4.6 The Heating System Analysis In Chapter II and III the mathematical model of the room was developed and validated using an experimental room model. In this chapter

TR T -T I Figure e.11. Analog Compuber Simulation for 2 he Simplified Room Model.

-88the remainder of the heating system mainly the heat exchanger, the thermostat, the air duct and the gas valve have been modeled mathematically. An analog computer setup was also constructed for each component. In the next sections the mathematical and computer models developed previously are utilized to study the interactions and dynamic behavior of these elements. A series of measurements were made with the computer model to determine the affect of thermostat hysteresis, thermostat time constant and furnace time constant on the cycling amplitude and period. The cycling amplitude is defined here as the difference between the maximum and minimum room temperatures when the system is operating. Similarly the cycling period is the time necessary for the room temperature to return to an original value while the heating system is in operation. For this definition, the system is assumed to be under constant outside environmental conditions. The tests and their results are described in the following sections. 4.6.1 Effect of Thermostat Hysteresis on Cycling Amplitude and Cycling Period In analyzing the effect of thermostat hysterisis, the outside temperature To is set at 20~F, and the furnace heat output is adjusted to give room temperature swings around 70~F. A sudden change in outside temperature is then applied. Simulating these conditions on the analog computer model for the heating system, recordings of the room temperature were obtained for a set of thermostat time constants and a set of furnace time constants. Samples of these recordings are shown in Figures k.12 and 4.13. From the information gained by these samples, Figures 4.14 through 4.17 were obtained.

-89Time Scale: 4 min/mm Room Temperature Scale: l~F/cm,,4 4 70 \.. a T 1 mm' 1.6 m q = Q... I T A O t 7 o l I i ii 1i,!.~.~, 7~ F:!!! i- i~ j:! j j'j>..il,. i. L L LLl time T F T - 1 minm, T = 1.6 minm q = 0.50F Figure 4.12. A Sample of Recordings. r.: ii:.i:.:- t.... time TT = 1 min. TF = 1.6 min. q = 0.3~:F -~ 70op!t iOF! 1, ii, 0 liil ~ ftii' i! I''A I-1.. 7' -t:~ i!l!ii i!F i i rti ~:i i I L "T = 1 min,'F ='1.6 min, q = 0.3.F Figure 4.15. A Sample of Recordings.

-90~ 1.6 ~~,,~, 0 z I 0.4 0 0.1 0.2 0.3 0.4 0.5 HYSTERESIS q Figure 4.14. Cnventional System's Cycling Amplitude-Hysteresis Curves with TT as a Parameter. 40.8 0.6 2 0.0 40 12..... ~~ 0 01 0.2.3 0.4.5 CC 24'T MN.. 20HYSTERESIS, U 16 12 Figure 4.15, Conventional System's Cycling-Period-Hysteresis Curves with TT as a Parameter.

-911.6 1.4 LL 0 1.2 UJ 1. 0 r Xy 0.8 _ _,IFr z 0.6 O 0.4 0.2 0 0.1 0.2 0.3 0.4 HYSTERESIS oF Figure 4.16. Conventional System's Cycling Amplitude-Hysteresis Curves with TF as a Parameter. 40 36 32_____ LL 24 0 0 0 0 0.1 0.2 0.3 0.4 HYSTERESIS OF Figure 4.17. Conventional System's Cycling PeriodHysteresis Curves with TF as a Parameter.

-92Inspection of these curves reveal the cycling amplitude and period increase with an increase of hysteresis. 4.6.2 Effect of Thermostat Time Constant on Cycling Amplitude and Period In order to investigate the effect of thermostat time constant, the heat output of the furnace is adjusted to give temperature room swings around 70~F. The furnace time constant is held fixed. The outside temperature is at 20~F and then is given a sudden change. From computer results it is seen that as the thermostat time constant decreases, the cycling amplitude decreases, for all values of hysteresis. It is also observed that with an increase in thermostat time constant, the cycling period increases, for all values of hysteresis. These are illustrated in Figures 4.18 and 4.19. 4.6.3 Effect of Furnace Time Constant on Cycling Amplitude and Cycling Period In studying the effects of varying the heat exchanger time constant, the thermostat time constant is held fixed. The furnace output is adjusted to give room temperature swings around 70~F. The outside temperature is set ao 20~F initially, and then given a sudden change. By simulating the model developed and validated in Chapters II, and III of a typical heating system on an analog computer, room temperature was determined and the response was recorded. Samples of these recording are shown in Figures 4.20, 4.21 and 4.22. In all the following cases a particular value of the furnace time constant and hysteresis were choosen, and the resulting values of the cycling period and amplitude were observed. The furnace time constant took on values from 1 to 5

-931.8..~ La. 1.6 -Ief Z44. O:3 -j 1.2 " — _03 g 0 OF O.8 0.6 0.4 0 I 2 3 4 5 6 THERMOSTAT TIME CONSTANT, TT, MIN Figure 4.18. Conventional System's Cycling Amplitude -TT Curves with q as a Parameter. 36 bO 28 Of 0 0 0 I 2 3 4 5 6 THERMOSTAT TIME CONSTANT T TT, MIN. Figure 4.19. Conventional System's Cycling Period -TT Curves with q as a Parameter.

-94. llil-till:-l: Time Scale: 4 min/mm:!::: 1: Temperature Scale: ~F/cm I TF = 4 min q = 0, Tp = 1min Figure 4.20. A Sample of Recordings. ^~~~~I~~~~~~~~f~!ifi iiii',1!~i iil i;ii7;iil 7 Tp = 5 min q = 0, Tr = 1 min Figure 4.21. A Sample of Recordings. I::::~:::ll ^ i;". ~~~time~ I..............:, ""....' " TF = i min, q = 0, TT = 1 min Figure 4.22. A Sample of Recordings.

-95minutes and the hysteresis from 0 to 0.3~F. From these results the cycling amplitude and the cycling period were plotted as a function of the furnace time constant. These curves are shown in Figures 4.23, and 4.24. From Figure 4.23 with the exception of the zero hysteresis case, it is seen that the cycling amplitude increases as the furnace time constant increases. As the hysteresis and the furnace time constant take on larger values, the rate of change of the cycling amplitude increases. In the case of zero hysteresis, the cycling amplitude increases until obtaining a maximum value, and then starts to decrease. From Figure 4.24 it is seen that the cycling period and its rate of change increases as the furnace time constant increases for all values of hysteresis. 4.7 Evaluation of Results The studies described above are the first phase in the application of the microscopic mathematical modeling of the heating system to the determination of space heating characteristics and requirements. The importance of the separate influeces on the system are seen at once by utilizing the model. The test performed with the computer model which was described in previous sections reveals what sort of information the model can best be expected to give. From Section 4.6.1, it is concluded that the room temperature cycling period and amplitude decrease with a decrease in the thermostat hysteresis. In addition it is shown in Section 4.6.2 that for a small cycling amplitude a thermostat with a small time constant is required, From Section 4.6.3, it is apparent,that, for a given value of

-961.6 1.4 o 1.2 - ~ - ~ DI~~~~~~~~~~~~~ / < 0 0.4 >" 0.8..... 0.2 1 2 3 4 5 6 FURNACE TIME CONSTANT, TI, MIN. Figure 4.23. Conventional System's Cycling Amplitude -TF Curves with q as a Parameter. 36, 32 Z 28 24 0 a20 o ~I 0 0 I 2 3 4 5 6 FURNACE TIME CONSTANT, T', MIN. Figure 4.24. Conventional System's Cycling period -T' Curves with q as a Parameter. O 2 3 4 5 6~~~~~~~~~~~~~~~~~~~~~~~~~01

-97hysteresis, the smallest cycling amplitude is achieved by selecting a furnace with the lowest possible time constant, with the exception of zero, hysteresis. The zero hysteresis is of no practical importance because such a control is not physically realizable. If in a given system, the hysteresis is relatively small, the furnace time constant plays a minor part in determining the cyclic amplitude. As a consequence of selecting a small cycling amplitude, the cycling period is small. In practice, a very small cycling period is not recommended since this will cause excessive deterioration of the gas valve. Therefore, a compromise between the cycling period and amplitude should be achieved. The results obtained in this chapter by studying the effect of parameter changes on the heating system do agree to a great extent with many experimental observations in the industrial field of heating systems. For example, a designer always attempts to devise a thermostat with the lowest possible hysteresis and time constant. He always tries to design a thermostat with a large UA product, and a small product of me. Efforts in the industrial field are presently being expended to replace the bimetallic thermostat with a solid state substitute because of its small hysteresis and time constant, The results in this chapter indicate what modifications can be expected from a change in heat exchanger time constant with a various combinations of thermostat time constant and thermostat hysteresis. That is to say, the mathematical model derived for the heating system can provide information about the relative effect of the three system parameters on the system performance. It can,moreover give this data, when additional variables such as a change in outdoor doncitions, or addition of a new

-98thermostat is introduced. These results are also in general agreement with many of the industrial heating system's experimental observations.

CHAPTER V OPTIMIZATION OF THE HEATING SYSTEM In the preceding chapters an adequate model has been developed to represent the forced-air gas heating system. In this chapter, the mathematical model developed previously is used to formulate a method for optimally controlling the heating system in accordance with a prescribed performance criterion. The optimal problem is one of reducing the room temperature deviation from a prescribed reference value to zero, while at the same time minimzing the value of some predetermined performance or cost functional J. The development of the optimal control law proceeds in essentially four steps. The first step consists of a reformulation of the mathematical model in a form which is more suitable for the application of optimization techniques. In the second step, an optimization criterion is defined which incorporates the main objective of minimizing room temperature variations with respect to a a prescribed reference temperature. In the third step, the particular optimization technique best suited for the optimization problem is chosen. Finally, in the fourth step, this optimization technique is utilized to construct the optimal control system. 501 Reformulation of the Mathematical Model, It is convenient to begin by reformulating the equations representing the fixed portions of the domestic heating system in the state variable formulation; i.e., in terms of a set of first order differential equations. -99

-100First, consider the equations describing the dynamics of the simplified gas-fired furnace as given in Chapter IV by(4.15) and (4.16). These equations may be expressed as follows: O = a7Te + Ta + a8TR Te = alTe + aTR + a3Tf where al, a2, and a3 are constants of the particular furnace. The state variable diagram for the furnace is shown in Figure 4.6. It should be noted that a number of assumptions have been made which render these constants independent of temperature and therefore, independent of time. Of course, this leads to a time invariant system, For the simplified hot air duct it was shown in Chapter IV: Ti = lTa + 2T5 where ~l = (r-i)/(-rl) and ~2 = 2/(r+l) A state variable diagram corresponding to these equations is shown in Figure 4.8. We now turn to the modified room model presented in Chapter III which includes the domestic space and the walls, and is described by the following dynamic equations:

-101QPCpTi - (k+QPCp)TR + kTw = pCpVTR -CR(R+Ro)TR + C(Ri+R)Tw - (2R +Ro)TR + (2R+Ro+Ri)Tw - RiTo = 0 Rearranging the equations representing the fixed components (i.e., room, air duct and heat exchanger) of the heating system, results in: = a T +a T +a (51) R 11 R 12 w 13 e Tw = 21TR + a22Tw + a23Te + m2 (5.2) Te = a1TR + a3Te+ U ( ) where the aij are constants which may be expressed in terms of t1, 2 Q Q P, Cp, k, V, R, Ri, and C. From the assumptions made in Chapters II and IV it follows that these aij values are time invariant. In addition, the variables u3 and m2 appearing in Equations (5.2) and (5.3) above are defined by: U3 = a3T Ri m2 = TO C(Ri+R) The set of linear time invariant differential equations expressed in Equations (5.1) through (5.3) may be expressed in the following single vector differential equation: x = Ax + u + m (5.4)

-102where: TR. ~ - 0 x = Tw u m =, m m e-3 0 The vector u- in Equation (5.4) is called the controlled input, the vector m is called the uncontrolled input, and A is the square matrix which specifies the linear time invariant system. The matrix A is given by: 11 12 113 A a21 a22 a25 a51 a32 a33 3- 1 32 where a32 0 The vector Equation (5.4) represents the state variable formulation of the dynamic equations describing the domestic heating system. 5.2 Formulation of the Variational System In this section a variational system is formulated to represent the heating system as referred to some equilibrium position. First, assume that the controlled input u and the uncontrolled input m are such that the system is operating in an equilibrium condition, in other words x = 0. In this case any disturbance which occurs in the system, for example an opened door, entering people, additional lighting, etc., causes a deviation in x from its nominal or equilibrium value.

-103To obtain the equilibrium values, set x = 0, therefore: AxO + uO + m 0 (5.5) The zero subscript here refers to the equilibrium vectors. In order to maintain a desirable room temperature which is a component of the vector x, it is evident that u the controllable input to the system in the equilibrium state takes on some value u. To determine the value of uo required, consider Equation (5.5): 11 12 T3 B [ll a2 aL i T i a21 a22 a23 = o 0 m2 0 Two0 21 ~ at e U3L By appropriate manipulations this equation becomes: a12 a13 0 TWol ta1 0 a22 a23 = - a21 TRo - m (56) From Equation (5.6), it can be seen that knowledge of the desired room temperature TR and the outside temperature expressed by m, determines the equilibrium values of the controllable variable Tfo and the various state temperatures. Let x = X +x, k = X+ 5x and u = uo + bu where 5x 6u, iu and 5x represent deviations from the nominal values of x, u, and x respectively.

-104Since the equilibrium point satisfies; the matrix differential Equation (5.4): xo = AXo + U + m (5.7) Also for any deviation 6x, bu, and 6x: (xo + x) = A(xo + 6x) + uo + u + m (5.8) From Equations (5,7) and (5.8) the following matrix differential equation for the deviations can be expressed as: 5x = A6x + 6u Now let 5x = y bu = v then y = Ay + v (59) This latter equations is the conventional well-known linear first order matrix differential equation. The components of the vectors in this equation represent variations in the state of the system. It is to be noted that in this case this new variational system (5.9) is valid for large swings from equilibrium since the original dynamic system is linear, The optimization problem now reduces to one of starting from some initial disturbance yo and driving the system (5.9) to the equilibrium state and at the same time constraining the original system to operate in such a way so as to minimize the value of some cost functional J(y,v).

-1055.3 The Optimization Criterion For a given dynamic system, the optimization criterion may assume a variety of forms depending on the requirements to be met. This is an important step in the design of the heating system since it determines to a large degree the nature of the resulting optimum controller. The main objective in the optimization of a gas-fired forced-air heating system is to reduce and penalize room temperature variations due to distrubances, and is primarily used here to define the optimization criterion. It is traditional in approximation problems arising in many fields to use a quadratic criterion. The square penalizing here will discriminate heavily against occasional large room temperature variations. This criterion expresses the philosophy that the only important aspect of the control action is following the desired output temperature and only imperfections in following are to be penalized. This philosophy is well justified as long as the type of control used does not have any significant physical limitations. In the case of a gas-fired heating system, physical limitations are imposed by the heat exchanger.(34) There is a minimum amount of energy input that the heat exchanger can give due to the characteristics of the gas burner. (4) There is also a maximum amount of energy which a heat exchanger can give to the system. Therefore in order to avoid the saturation effect of the heat exchanger, a term in the square error criterion is added that is proportional to the square of the control signal. Having these two factors in mind, the optimization criterion for the forced-air heating system can be represented as follows:

-1063 T J(yv) = Z f [q2(a) + v2(a)]do (5.10) n=l o n n where J(y,v) is the error criterion to be minimized, a is a dummy time variable, T is the period over which the minimization takes place, and qn(a) is defined as: ql -- 1 0 ql(a) 1 0 O Yi(a) q(a) - q2(a) = 0 O Y2() q3(a) 0 0 0 V3() The first term in the integrand of the quadratic criterion (5.10) represents the penality on the room temperature variations, and the second term is introduced to avoid saturation of the control signal and also for power limitation. 5.4 The Optimal Control Law In recent years, many elegent mathematical techniques have been presented evolving the theory of optimal control. Most of them are based on the calculus of variations, the maximum principle and the dynamic programming. Excellent contributions to the theory of optimization can be found in the literature. Much of the work, however is of a mathematical nature, since it was written by mathematicians. In this section one of those optimization techniques mentioned above which is most suitable is used to solve for the optimal controller of the optimization criterion defined in section (5.3). In many of the cases treated by the calculus of variations and the maximum principle techniques, a control law of an open loop nature results.(23' 4) This is not desirable in the heating system case, since

-107unwanted distrubance signals at the output of the system (room temperature) will cause errors. Meanwhile, the method of dynamic programming applied to this linear time invariant heating system is guaranteed to provide a closed loop or feedback control law(l4'25) for a given set of heating system parameters, which satisfies the optimization criterion defined in section (5o3). The Euler-Lagrange equation resulting from the application of the calculus of variations is always unstable and poses some difficulties.(14) This provides the dynamic programming method with another advantage, since it does not provide such difficulty. Also, if the optimization process is carried over a semi-infinite interval (as T -oo), difficulties will arise in using the calculus of (14) variations, and the maximum principle methods. () However, as will be seen, the method of dynamic programming overcomes such difficulty. For these reasons, the method of dynamic programming is thought to be the most suitable method for the optimization of the heating system under the optimization criterion represented by (5.10). Bellman's Dynamic programming is basically an optimization process that proceeds backward in time; that is, the solution is computed over the last interval of the process and successive solutions are computed for the remaining intervals of decreasing time until the total solution is obtained for the entire process. The optimization problem now at hand is one of starting from some initial temperature disturbance yo, and driving the system = Ay + v to the equilibrium state while constraining the original system to perform in such a way as to minimize the value of the cost functional J(y,v). Here the period of optimization is allowed to be

-108very large (i.e., T -oo), since the heating system has to be optimized over a long period of time. In order to apply the functional equation technique of dynamic programming, this optimization problem is embedded within the wider problem of minimizing: T 2 [2) + 2 I [() 4- Vn(a)]da n=l t subject to the heating system Equation (5.9) and the initial condition y(to) = yo, with t ranging over the interval (to,T). Let the minimum of this cost functional be: 3 T E(y,t) min 7, [q2(a) + v ]da (5.11) v n=l t n n Invoking the principle of optimality to Equation (5,11), the functional equation becomes: 3 t+~ 2 2n E(y,t) min { [qn() (a)]da+ E(y+ye, t+e) } v n=l t ( (5,12) which is reduced to the following expression (by integration and Taylor series expansion): E(yt) min v{ [qn(a)+Vn(a)] + E(y,t) + Z y E + e} +A (e) v n=l n n n=l Simplifying: 5 3 2 3 mn {l[q(a)+vn(a)] + 1 Yn + + (E) = where Z (E) -40 as e - O

-109Therefore 2 2 bE min { [q2(a) + vn(a)] + - } = 0 (5.13) v n=l n=l nYn t The minimizing control signal vector v*(a) is obtained by minimizing the sum of terms within brackets with respect to each signal of the control vector. Minimizing now with respect to v3(a), the only non zero component of the vector v(a), and keeping in mind the relation between the vectors q and y, therefore: 2v* + E = 0 3 aY3 where v* = optimum control signal. Consequently, the condition for minimum error is: v = 1 E (5.14) 3 2 ay In order to determine the optimum signal v*, -E for minimum error 3 ay3 should be determined first. By substituting Equation (5.,14) and the value of q in terms of y into the functional Equation (5.13), the condition for minimum error becomes: 2 + (E )+ Z + o (5.15) I 4 y n=l = 4Yn t As seen from Equation (5.15), the condition for minimum erroris in a partial differential form. To solve such an equation a power series solution is assumed, then the coefficients in the series are found by direct substitution.

-110Since the integrand of the error criterion function is a quadratic expression and the dynamic system is linear, the minimum error function E(y,t) is also quadratic and can be written as(33) 3 3 3 E(y,t) = k(t) + Z km(t)Ym(t) + Z Z kmk(t)Ym(t)Yk(t) m=l m=l k=l (5.16) where kmn(t) = knm(t) where k(t), km(t), kmn(t) are the parameters to be determined from Equations (5,15) and (5.16). By partial differentiations of Equation (5.16), (aE(y,t))/(5yn) and (aE(y,t))/(6t) are written as follows: aE(y,t) = k (t) + 2 k (t)y(t) (517) 6yn n m=l n M and k'(t) + E k~(t)Ym(t) + E E kmk(t)Ym(t)Yk(t) =6t m=l m=l k=l (5.18) If these partial derivatives are substituted into Equation (5.15) the condition for minimum error becomes: 2[3Zky + Z knmY]2 3 y+ [k +2 k my + k + E k +m k' IY 35 m=l m=l mm m=l k=l mkYmYk (5.19) 3 3 + [k y + 2y E k y ] 0 n=l n n n m= nmm The condition for minimum error expressed by (5.19) is satisfied for all finite values of yn(t), assuming the k-parameters are independent of yn(t), only if each of the coefficients of the constant term, yn(t), and yn(t)Ym(t) in Equation (5.19) vanishes, where n,m = 1,2,3. Therefore by equating the coefficients of the constant

-111term, yn and YnYm each equal to zero, the following simultaneous first order differential equations in the k-parameters result. f (k'kl k2 k3, kll k22, k33, klk2kl3, k23) = k' f5(kklk2,k3kllk22,k335k12,k1,k2 -) = (5.20) flo(kklk2,k3,kll k22,k33,k12, k13,k23) 23 where: fl, f2.... fl are in general nonlinear functions of the k-parameters. The primed k's refer to the derivatives of the k-parameters with respect to time. This method of assuming a solution leads to the reduction of the problem of solving a partial differential equation to the problem of solving a set of first order ordinary differential equations. The boundary condition for the k-parameters are deduced directly from the required boundary condition on the minimun error function. From the expression for minimunm error function for t = T, the boundary condition is E(y(T), T) = 0 which means that k(T) = kn(T) = knm(T) = 0 (5.21) The problem becomes now one of finding the optimum control system of a one-point boundary value problem. The parameters of the

-112optimum control system k(t), kn(t), kmn(t) where mn = 1,2,3 can be determined from the set of ten differential Equations (5.20) with boundary conditions given by (5.21). It is to be noted that the number of parameters are ten and the number of initial conditions expressed by (5.21) are ten. The solution of the set of differential Equations (5.20) as T tends to oo, must assume steady state. If the k-parameters assume steady state values, then the differential equations given by (5.20) reduces to a set of algebraic equations. Therefore, when the dynamic system is time invariant, the error function is quadratic, and the optimization process is carried over a semi-infinite time interval, the parameters of the optimal qontrol law become time-invariant. Since the heating system is to be optimized over a semi-infinite time interval for a quadratic optimization criterion, the minimum value of the error functional becomes only a function of the state variables and not a function of time and expression (5.11) becomes: E(y) = mn l r[q(a) + v (u)]dJ (5.22) v n= t' n n Invoking now the principle of optimality to Equation (5.22), the functional equation becomes: E(y) m in ~m t+c 2 2 E(y) { [qn()+n()]d E(y+(y+e)} (5.23) _n=l t which is reduced to the following expression by integration and Taylor series expansion: - E(y) { [q2()+v2()] + E(y) + E Y n _e + A(E) v =1 n n nyn

-113Since A (e) -,O as e -O, therefore: 3 2 23 min [q (a) + n()] + n } = (5.24) [ n- ( c) +vcf I + F''n n v n=l nn=l Applying now the minimization process to Equation (5.24) similar to the minimization of the functional in Equation (5.13), and having in mind that (3E)/(3t)= 0 for the semi-infinite time interval case, therefore: = _ (5.25) 3 2 ay3 which defines the optimum control vector, and the condition for minimum error (5.15) becomes: 2 1 aE )2 aE ( ^,L+ + y Z 10 (5.26) 14 + n =i n y Expression (5.16) for the semi-infinite time interval case takes the following form: 3 3 3 E(y) = k+ Z kmy(t) + Z Z kmkym(t)Yk(t) (5.27) m=l m=l k=l where k, km, and kmn where m,n = 1,2,3 are fixed constants. Differentiating Equation (5.27) now partially with respect to y, n = 1,2,3 gives: kn + 2 kny(t) (5.28) aYn m=l and E, = k~ + 2[k1y1 + k2Y2 + k3y] (.29) ay,~~~~~~~~~~~(.9

-114By substituting (aE)/(ay3) from (5.29) into (5.25) gives: k = _ 3- k y k y (5530) 53 2 51 522 3Y3 Therefore it is necessary to determine the parameters k3, k31 ^ k32 and k33 to define the optimum control signal. Substituting now from Equations (5.28) and (5.29) into the condition for minimum error (5.26), and also using the vector matrix differential equation y = Ay + v, the following is obtained: 3 3 3 3 3 y 1 [k2 +4k3 3ym+4( m m) a mnm) [kn+21ksnys] 0 1 4 m=l m=l n=l m=l (5.31) By equating the constant term in Equation (5.31) to zero, since this equation is satisfied for all values of yn(t), the following is obtained: k3 - 0 (5.32) Similarly for the coefficient of y 3 -k3kml + amnkn 0 m 1,2,3 and since this is true for all finite values of a, therefore mn kl = k2 = k3 = 0 For the coefficient of y: 2^ J 3 n (3 1 - kl3 + 4 nZl an k = 0 (5.33) 15 ~n~ll

-115For the coefficient of y2: 2 3 - k23 + 4 a2n k2n = (5.34) 2 n=l For the coefficient of 2: 3 - k2 + 4 Z a k = 0 (5.35) 33 n=l 3n 3n For the coefficient of Y1Y2: 3 - k3 k2 + n (ln kn + a kin) 0 (536) For the coefficient of Y1Y3: ki k33 + nZ (aln kin + 0n kin) 0 (5.37) 13 23 n=l in 2n 2an k1n For the coefficient of Y2y3: 3 2k3 k33 + n (a n k3n a3n kn) (538) For the coefficient of y2y3 ~ -k k + Z (a k+ a kn -- 0 (5 8) 25 55 ni 2n= 2n 5n 2n Equations (5.33) to (5.38) are in general nonlinear algebraic equations the parameters kmn ) mn = 1,2,3 and requires a digital computer for solution. In the next chapter a solution of these parameters for a particular heating system on the digital computer will be shown. Equation (5.30) now becomes: V*= -[k y +k y +k y] (39) 3 5311 32Y2 33 3 The control function v*, derived here is referred to as the optimum control law,

This optimal control scheme for the variational system may be combined with the equilibrium system developed in Section (5.2) to obtain an optimal feedback system for the heating process. In block diagram form the system may be schematically represented as shown in Figure (5ol). In this diagram a controller is provided which compares the values of the environmental state and the desired state and commands the appropriate equilibrium input. It should be noted that the number of feedback loops is equal to the order of the heating system; it is noted also that the feedback signals are the measurable state variables. Thus, an optimal heating system for a defined quadratic cost function has been developed which has the desirable property of providing feedback loops to account for disturbances in the system. This optimum law will be applied to a particular heating system in the next chapter to develop an optimum controller. The optimal heating system is then simulated on an anlog computer to study its behavior.

5 0 U x -t +~~~~~~~~~~~~~~~~~~~ X =Ax+ u + m SENSORS H F i g u r e 5. 1. B l o c kmDi a g r a m-of' t h e O p t i m u m H e a t i n g S y s t e m. Yl k1 Y3 k33 Figure 5.1. Block Diagram of the,0ptimum Heating System.

CHAPTER VI THE ANALOG COMPUTER SIMULATION OF THE OPTIMAL HEATING SYSTEM In Chapter V, an optimum controller based on a quadratic error criterion was defined for the heating system mathematical model developed in this study. For optimization over a semi-infinite time interval, the feedback loops of the optimal heating system were shown to be timeinvariant, and consisting of time invariant gains multiplying the measured state-signals. However, the parameters that define the optimum controller must be determined for each particular heating system by solving the set of Equations (5.36) to (5.41) using either an analog or a digital computer. In this chapter, the particular heating system components (including: the room, the air duct, and the heat exchanger) considered in Chapter IV are combined to represent the fixed portion of the dynamic system. With this as the given system, the optimum controller is then determined according to a quadratic error criterion with the results developed in Chapter V. The parameters of the optimum controller are determined through the use of the 7090 digital computer facility at the University of Michigan. The resulting optimum heating system is then simulated on the electronic analog computer in order to investigate its performance, 6,1 The System Equations: The following equations were obtained for the heating system components (heat exchanger, air duct, and the room) studied in Chapter IVo TR = - 0.191 TR + 0.0422 Tw + 0.097 Te (6.1) -118

-119Tw = 0.2278 TR - 0.0974 Tw - 0.097 Te + 0.0184 To (6.2) Te = 0.25 TR - 0.489 Te + 0.239 Tf (6.3) These linear, time invariant, first order differential equations represent the dynamic process whose performance we wish to optimize. Equation (6.1), (6.2), and (6.3) may be written in the vector form to yield: x = Ax + u+ m, where 1al a12 a13 -0.191 0.0422 0.097 A = a i21 a22 a23 = 0.2278 -0.0974 -0.097 31 32 3a3 { j.25 0 -o.489 u= 0 Tf, m = 0.018 To 0.239 0 The variational vector matrix differential equation as derived in Chapter V for this heating system now becomes: y = Ay+v where v1 0 v = v = Tf V3;.0.239 where v1 = v2 = 0, and v3 = 0.239 6Tf

-120and Y1., 6TR Y = Y2 =- where Yl= 6TR Y - T y2 y e 6.2 The Parameters k31, k32, and k3 of the Optimum Controller The square matrix A defined in Section 6.1 determines the system under consideration, and therefore the k parameters of the system as defined by Equations (5.36) to (5 41) may be written as follows: 1 k2 0.388 k -0764 k + 0 1688 k - 0 (6 4) 31 51 11 12 - k 02 - 0.388 k32 + 1.112 k12 - 0.3896 k22 = 0 (6.5) k33-13 1 (6 6 1- k132 - 0.2884 k12 + 0.0422 k22 + 0.097 k32 + 0.2278 kll - 0097 k1 = 0 (6.7) - k31k33 - 0.68 k31+ 0.0422 k2 + 0.097 k33 + 0.25 kll = 0 (6.8) - k2k33 + 0.2278 k31-0.5864 k32 - 0.097 k33 + 0.25 k12 = 0 (6.9)

-121The University of Michigan Control System Algorithm Program employing a 7090 digital computer was used to solve for the k-parameters. This program was basically obtained from IBM, with some modifications added. The modified program is entitled CSAP and is currently available at the University of Michigan Computing Center Library. This program appears as a subroutine on the system disc and may be entered simply by calling CSAP. Once the program has been called, it will function exactly as described in the user's manual. () The solution for this particular system is: k = 1.4449, k22 = 1.3273, k33 = 0.2065 k2 = 0.805, k1 0.4467, k = 6.36. From Equation (5.39) the optimal control signal becomes: v* = - (0.4467 yl + 0.36 Y2 + 0.2065 y), or 0.239 8Tf = - (0.4467 sTR + 0.36 sTw + 0.2065 bTe) (6.10) Hence, the block diagram of the optimum heating system using this control law follows as shown in Figure 6.1. 6.3 Analog Computer Simulation of the Optimum System: 6.3.1 Fixed Part of the Heating System By making the time scale such that every second on the electronic analog computer corresponds to four minutes of real time, the differential Equations (6.1), (6.2), and (6.3) representing the optimal heating system become:

xH.23^ t ^.239 I9 I I~1 ^^ ^""^ x ^^ ^^^^ v~~x ~x~ - Ax+ u + m — ) Sensors.4467 r').36 8T, I~~~~~~~~~~~~~~~~~~~~~ Figure~ 61 BlcDarmoteO.20e6 Figure 6.1. Block Diagram of the Optimum Heating System.

-123T = - 0.764 TR + 0.1688 Tw + 0.388 Te Tw = 0.9212 TR - 0.3896 Tw - 0.388 Te + 0.0736 To Te = TR -1.956 Te + 0.956 Tf The temperature scale is selected in such a way that one volt is equivalent to 20F. As a matter of fact, the model developed previously in Chapter IV is used here, 6.3.2 The Simulation of the Optimum Controller In Section 6.2, the optimum control law was shown to be: 0.239 sTf = - (0.4467 5TR + 0.36 5Tw + 0.2065 STe) where 6Tf, TTR, 8T, and 6T are the variations of temperatures from equilibrium values, and are defined as follows: ST T T 5TR T TRo ) 8Tf = T f Tf TR TR- Ro 6Tw = Tw - Two' Te = Te - T Therefore, to generate 8Tf, the variational control signal, it is necessary to first generate the equilibrium values of the temperatures, Tfo TR Two and Te 6.3.3 Equilibrium Values: It was indicated in Chapter V that given a set room temperature TRo, and an outside temperature To, the equilibrium values of the control signal Tf and the state signals Two and Teo can be generated by the following procedure:

-124For equilibirum conditions: TR Tw Te = 0 and therefore under equilibrium conditions (6.1), (6.2), and (6.3) become: T = - 1.045 TR + 2.05 Te (6.11) Te = - 0 43 T + 1. 97 TR (6.12) T 237 TR - T + 0,19 To (61i3) wo Ro eo ~ wo 0 0 From these latter equations it follows that having TRo, and To as set inputs, the equilibirum values Tfo Teo, and Two may be generated on the analog computer as shown in Figure 6.2. Having established the equilibirum values generated as shown in Figure 6,2, they may be now combined with the fixed portion of the heating system and the optimum controller shown in Figure 6.3 to provide the analog computer simulation for the optimum control system, 6,4 Analysis of the Optimum Heating System: In this section, three runs on the optimum heating system simulation are considered. These runs were made for the purpose of analysis and present the results of selecting a quadratic performance criterion, The situations considered in the three runs were: 1o First Case: The room temperature TR is set at 70~F, and the outside temperature is initially set at 20~F. The system is therefore initially in the equilibrium state of:

-TO TWO C I~ 0..TR 197 10 *ec A 0 ^o ^1 T o Figure 6.2. Analog Computer Simulation for the Generation of the Equilibrium Values.

_Te R TTR /T ~~0 3~~~ ~0 -Te ^1 -^w ^'TRo I ~~~~____ IH Fgr65 mogopt DigmfrhOpmmHtnSse-mTw Two Figure 6.3. Analog Computer Diagram for the Optimum Heating System,

-127rTRO 70 x~o= Two = 52,9 F, Teoi 115 Tf = 163.80F The outside temperature To is then suddenly changed from 20~F to 0~Fo For these conditions the room temperature TR, the surface wall temperature Tw, the heat exchanger wall temperature Te, and the control signal temperature Tf were recorded. These temperature responses are shown in Figures 6.4, 6,5, 6.6, and 6.7 respectively, From Figure 6.4, it is evident that the room temperature TR decreases gradually from the time the disturbance occurs until the time when the variation 8TR becomes -0.15oF; a total of 16 minutes, After this, it begins to increase at a slower rate back toward its original value. In 40 minutes, the room temperature attains the value of 69 90~F This is expected, since the optimization criterion was considered over a semi-infinite time interval, The optimum controlsignal Tf as seen in Figure 6.7 increases gradually from the time of the drop in the outside temperature. This effect occurs to compensate for the heat loss caused by'sudden distrubanceo In Figure 6,5, it is noted that the surface temperature of the wall initially falls rapidly to.45.37F, then it gradually begins increasing until it reaches 49~F, Figure 6.6 indicates the effect of the disturbance on the heat exchanger t.emperature Te, This temperature initially drops to about 109~F because of both the decrease in room temperature and the decrease in surface wall temperatureo It then begins to gradually increase until it reaches within 2,4~F of its original value. This is caused by the increase in the flame temperature.

-128In the Graphs below: All Variations are Measured from the Original Equilibrium Position. 2 OF TR.1 0 -.2 Figure 6.4. Room Temperature Variation Response for Case 1. 10 5Tw OF 5 0.10 Figure 6.5. Surface Wall Temperature Variation Response for Case 1 10 -10 Figure 6.6. Heat Exchanger Temperature Variation Response for Case 1. 20 - Tf oF 10 -10 -20 TIME IN MINUTES...! I,I.I I I I I I I I 0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 Figure 6.7. Flame Temperature Variation Response for Case 1.

-1292, Second Case: For this case the room temperature is set at 70OF and the outside temperature is initially at O~F. The equilibrium values are: Roro x - Twol = 49 OF U 117 [eoJ Tf = 182,80F o The outside temperature then rises suddenly to 40~F. Figures 6.8 to 6.11 show the state variables and control signal responses to this disturbance, The room temperature response is shown in Figure 6.8. This temperature increases gradually to 70.15~F then falls to 70.05~F. The wall surface temperature shown in Figure 6.9 rises as a result of the disturbance and then it decreases until it reaches 63'F, The heat exchanger temperature also rises by 60F and then will decrease gradually to 110F This result is illustrated in Figure 6,10o Figure 6,ll shows the optimal control signal, It is apparent that the flame temperature changes gradually to 153,5~F, which implies that the disturbance caused by an outside temperature rise from O~F to 40OF, decreases the flame temperature by 29,03F to keep the room temperature to within OolF of 700F. 3. Third Case: The room temperature is set at 70~F, and the outside temperature initially maintained at O0F, This makes the equilibrium values of the state variable:

-130In the Graphs below: All Variations are Measured from the Original Equilibrium Position..2 0 FR Figure 6.8. Room Temperature Variation Response for Case 2. -.2 OF 20 0 -20 Figure 6.9. Surface Wall Temperature Variation Response for Case 2. -40 _ 101- Te OF -10 Figure 6.10. Heat Exchanger Temperature Variation Response for Case 2. 30 oF Tf 15 C -30 TBME IN MINUTES.. I. I I I I. I. 1 1.I I.. I _, 0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 Figure 6.11. Flame Temperature Variation Response for Case 2.

-131 TR 7O xo w 491 oF, Tf = 182.8~F o Now the outside temperature To is suddenly increased to 20~F. Figures 6o12 through 6,15, indicate the state variables, and control signal responses to such a disturbance, The maximum room temperature variation observed, 6TR, is 0.150F, and the final value is less than 0.10F. The control signal variation is seen to be -19~F meaning that the heat exchanger must supply a flame temperature Tf which is 19F less than the original value to compensate for a 200F outside temperature rise from 0~F. 6.5 Conclusion The results of this chapter indicate that the design specifications for the heating system can be met, The electronic analog computer simulation of the optimum heating system shows that successful automatic feedback heating systems can be achieved utilizing the quadratic error criterion. For the various outside environment conditions and disturbances, the optimum system appears to drive the room temperature TR to within 0,10F of the equilibrium or set value. The control signal Tf is driven to a value in each of the conditions studied which agrees with practical requirements as far as the level of the signal is concerned, In other words, the control signal does not assume upper values which cannot be implemented practically, and also does not take on lower values which a physical heat exchanger is unable to supply.

-132In the Graphs below: All Variations are Measured from the Original Equilibrium Position..2 oF | 6TR.1 0 ~ -.1 -.2 Figure 6.12. Room Temperature Variation Response for Case 3. 10 ~F 8Tw 5 -5 -10 ~ Figure 6.13. Surface Wall Temperature Variation Response for Case 3. 10 8Te 5 OF -0 -10 Figure 6.14. Heat Exchanger Temperature Variation Response for Case 3. 20 8T 10 OF -10 -20 TIME IN MINUTES I I I I I I.- I I I.,, I I I 0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 Figure 6.15. Flame Temperature Variation Response for Case 3.

-133If this system is to be implemented in practice, some practical problems have to be solved. The implementation of this system requires knowledge of the various system constants appearing in the matrix A and elsewhere. Some of these parameters may be difficult to measure however, an estimate can be reached by the knowledge of the construction details of the system. The other problem might come in adding complex equipment to measure the state variable Te which is the heat exchanger temperature. However, the remainder of the state variables will be easy to measure. Another practical problem arises from the fact that this system requires a furnace which is able to be modulated in a continuous fashion. Fortunately great effort is being expended at the present time in industry in order to produce such systems. A modulation system has been designed and tested at the Maxitrol Company in Detroit. Its superiority over an on-off system was demonstrated when tested in a model house, however, for economical reasons it cannot be used for domestic heating at least at the present time. However, disregarding for a moment the practical problems which arise, the system above does represent an optimum for the configuration and cost function chosen from the theoretical point of view. Thus it represents a bound or a standard with which conventional or sub-optimal systems may be compared. Also, for some specialized installations with rigid performance standards it may be feasible to utilize a system such as this. From an overall viewpoint, although the particular optimal system developed here may not be readily applied to current domestic heating systems mainly because of economical reasons, it seems that

-134valuable information might be gained by investigating sub-optimal adaptations of this type of system. In the next chapter a sub-optimal controller for the conventional forced-air domestic heating system is proposed based on results from the optimal control law derived in this chapter. This sub-optimal controller is simulated with the conventional heating system on the analog computer to compare results with the conventional controller, and the optimum controller,

CHAPTER VII A SUB-OPTIMAL SYSTEM In the preceding chapters the theory of optimal control has been used to evaluate the form of the control signal and the structure of the heating system for a quadratic error function. It was apparent from the complexity of the equations resulting from the optimalfeedback control system that a digital computer was needed to solve the equations. Therefore, before one attempts to obtain optimal solutions numerically one should first consider the possibility that a conventional controller might do almost as well as the optimal one. In addition, if the optimum system proves to be uneconomical, an investigation should also be made for the use of a sub-optimal controller. In a conventional heating system, the thermostat (system controller) senses the room temperature TR and sends a signal to the furnace which turns it on if the room temperature is lower than a set value, and turns it off if the room temperature is higher than this set value. The optimum control system derived in the previous chapter calls for sensing not only the room temperature TR but also the wall surface temperature Tw and the heat exchanger wall temperature Te These three signals are fed back to the heat exchanger through a set of feedback gains and thus continuously vary the output of the furnace, A fourth sensor is also necessary for monitoring the outside environment temperature To which determines the equilibrium operating values. Thus, to implement the optimum system in practice, four sensors and their associated feedback gains are needed to sense the state variables, -135

-136and continuous modulation must be added to the heat exchanger. In addition summing equipment is needed to produce the variations 6TR, 6Tw, 5Te~ and 6Tf. Thus the optimum system is in general expensive unless it is used for a specialized installation with rigid performance standards which are independent of cost. For domestic heating it would be particularly costly due to the additional expenses attributable to a continuously modulated furnace, four temperature sensors, and the associated circuitry making the optimum system impractical. Therefore, it is logical to attempt to devise a sub-optimum heating system, which (1) can be implemented practically, (2) is acceptable economically, and (3) yields improved performance compared to the conventional system. In this chapter such a sub-optimal control system is proposed. It is based on the properties of the optimal system developed in the previous chapter. This sub-optimum system is simulated on the electronic analog computer, and compared with a conventional control system and the optimum system. The chapter is concluded by describing a method for implementing the sub-optimal controller. 7.1 The Proposed Sub-Optimum Heating System The principal properties of the optimum heating system developed in Chapter 5 may be summarized as follows: (a) Three state signals, the room temperature TR, the surface wall temperature Tw, and the heat exchanger temperature Te are combined to form the control signal. (b) Three feedback gains are required to provide the proper weighting of the state signals TR, T, and Te

-137(c) An external sensor is necessary to sense the environmental temperature To, and generate the equilibrium values. (d) Summing equipment is employed to generate the variations 5TR' 5Tw, 5Te, and 5Tf (e) A heat exchanger with the capacity for continuously modulating flame temperature is required. For a conventional heating system the main properties can be summarized as follows: (a) A single signal, the room temperature TR, is used for control purposes, (b) A heat exchanger with only an on-off capacity is used. The characteristics of any practical sub-optimum heating system which are proposed should draw heavily from the properties of the optimum heating system. To develop a sub-optimum system incorporating the properties of the optimum and conventional systems, the following reasonings are used: (i) In present day forced-air gas heating systems a heat exchanger with continuous modulation equipment is considerably more expensive than one without such equipment. The improvement in the system due to adding such equipment versus the cost of the equipment is not acceptable at the present time to the consumer. Therefore it is essential that the sub-optimum system utilize an on-off heat exchanger. (ii) In both conventional and optimum heating systems room temperature TR is sensed, and controlled. Therefore, it follows that

-138room temperature TR should be controlled in any sub-optimum system which is proposed. (iii) In the optimum heating system, the wall temperature Tw is sensed and utilized for control purposes. This means that the Tw sensor attempts to predict environmental phenomena, by sensing the wall surface temperature, and uses this prediction to send a signal to the furnace to react before significant changes occur in the room temperature. Intuitively this principle can be applied in a sub-optimum system. Therefore, it is suggested here that the sub-optimum system include a sensor to monitor Tw and feedback a signal to improve variations in the room temperature. Tw the equilibrium temperature of the wall surface is set at a value corresponding to a room temperature of 70~F and the average outside temperature under consideration. (iv) The heat exchanger temperature Te sensed in the optimum system should bot be sensed in the sub-optimum for economical reasons. For the reasons stated above the proposed sub-optimum system has the following main characteristics. (a) Heat is applied from a heat exchanger in an on-off fashion. (b) The room temperature TR is sensed and fed back for comparison with a reference signal. (c) The surface wall temperature Tw is sensed and fed back with an appropriate gain factor. This proposed sub-optimum heating system is shown schematically in Figure 7.1.

x TR T Fgr71.l ok DarmShmtco h u-piu etn ytm Figure 7.1. Block Diagram Schematic of the Sub-Optimum Heating System.

-1407.2 Analysis of the Sub-Optimum Heating System The sub-optimum system described in Section 7.1 was simulated on the electronic analog computer for the purpose of investigating and analyzing its performance. The fixed portions of the heating system used in Chapter 4 were used in order to compare the different controllers developed thus far, namely; conventional, optimum and sub-optimum controllers. The analog computer simulation for the sub-optimum heating system is essentially the same as the one used for the conventional heating system developed in Chapter 4, with the addition of a feedback loop from Tw. Thus it is not presented here. In this investigation, the heat output of the furnace is adjusted so that the temperature of the air circulating in the heating system is 120~F when it is leaving the furnace during the on-period. During the off-period the temperature of the air is considered to be 70~F. The outside temperature To is held fixed at 20~F, and then allowed to drop suddenly to zero. Computer runs were made for the system with room air temperature TR feedback alone, and then with both room air temperature TR and surface wall temperature Tw feedback. Let kw be the feedback gain for Tw. In each run, the room air temperature response cycling amplitude and period were recorded. Computer runs were made for different Tw feedback gains. For each of these runs, the thermostat time constant, heat exchanger time constant, and the thermostat hysteresis were varied and plots of the cycling amplitude and period versus the feedback gain of Tw were obtained for various sets of parameters. The plots and samples of recordings of these investigations are shown in Figures 7.2 through 7.9.

-1417 TT = 0 4 min T = 1.6 min.5 -.4 q 0.18~F 3 q = 0.12~F.2 q = OOF.1, I,... I I I 0.2.4.6.8 1 1.2 Figure 7.2. Cycling Amplitude-kw Plots with q as a Parameter. 36 - T 0.4 min TF = 1.6 min 32 F q = 0.18~F 28- 24 820 q = 0.10~F 12, I.I.I. I w I. 0.2.4.6.8 1 1.2 Figure 7.3. Cycling Period-kw Plots with q as a Parameter.

-1420p In Figure 7.2 three curves are displayed. These curves indicate the dependence of the cycling amplitude on the feedback gain of Tw for three values of thermostat hysteresis, q, and for a thermostat time constant, TT, of 0.4 minutes and furnace time constant, TF, of 1.6 minutes. An examination of the curves shows that for the range of feedback gain considered, the cycling amplitude decreases as a function of feedback gain and then begins to either increase slightly or become constant. The case of a zero feedback gain for Tw represents the conventional controller. For a hysteresis value of 0.180F the cycling amplitude decreased from 0.630F to 0.330F when the controller was changed from a conventional form to a sub-optimal form with a unity feedback factor. This indicates an improvement of about 50% in the room temperature deviation. Figure 7.3 displays the corresponding curves for the cycling period. It reveals that by starting with a higher value of thermostat hysteresis, the cycling period decreases with feedback gain, and the rate of decrease of this period drops with the thermostat hysteresis. For the case of a zero hysteresis the cycling amplitude tends to increase slightly after it reaches a minimum. Figure 7.4 displays three curves showing the cycling amplitude versus the feedback gain for different values of heat exchanger time constant. The thermostat time constant was kept at 0.4 minutes and its hysteresis at 0.120F. Examination of the curves indicates that for heat exchangers of time constant 5 minutes, 2.5 minutes, and 1.6 minutes, the cycling amplitude is a decreasing function of the feedback gain of Tw. For values of feedback gain kw > 0.9 the cycling amplitude begins to be constant. Therefore, when a heat exchanger of a time constant of 5 minutes is used with the heating system and a conventional

-143- T = 0.4 min, q = 0,12F.6'r~-F = 5 min I'5 TF = 2.5 min.4 o Tp = 7F = 1.6 min.3 0.2.4 6.8 1 1.2 Figure 7.4. Cycling Amplitude-kw Plots with TF as a Parameter. T = 0.4 min, q = 0 12 OF 40 36 32 TF = 5 min 28 24 20 i ~~20 \ TF = 2.5 min 16a 12 48 I l I I I 0.2.4.6.8 1 1.2 ] Figure 7.5. Cycling Period-kw Plots with rF as a Parameter.

-144TF = 5 min 0.9 q =.12~F rr 0.8 " $ -- 15 -= 0. 9 min + 0.7,~ ^^* ^ TT -= 0. min ^ 0~.5 ^^T" "T.-___ = 0.3 min 0.4 0.3 0 0.2 0.4 0.6 0.8 1 1.2 Figure 7.6. Cycling Amplitude-kw Plots with TT as a Parameter. 40 TF = 5 min F 36 q =.12~F 1F 32 3 2 ^ 24 4^ 7=0,4 min. 20 0 0.2 0.4 0.6 0.8 1 1.6 Figure 7,7. Cycling Period-kw Plots with TT as a Parameter.

-145~. -^..~... Temperature Scale: O. ~F/cm *g: v* *i.. Time Scale: 1 min/mm E- fI^^!_t ifi -time T = 04 min, T = 1.6 min q = OF = 0.8 Figure 7.8. A Sample of Recordings. TT = o.4 min, TF = 1.6 min, q = 0~F k=0 Figure 7.9. A Sample of Recordings.

-146controller, it shows a 0.657~F cycling amplitude and 0.55~F cycling amplitude when used with the suggested sub-optimal controller with 0.725 gain, i.e. it shows an improvement of approximately 19% when the sub-optimal controller is used. Figure 7.6 shows the corresponding cycling period set of curves. At lower values of heat exchanger time constant, the cycling period decreases gradually, and as the heat exchanger time constant increases, the room temperature cycling period decreases and subsequently begins to increase at a very small rate. Figure 7.6 displays a graph indicating the cycling amplitude versus the feedback gain, when the thermostat time constant is varied as a parameter. The heat exchanger time constant was held constant at 5 minutes, and the thermostat hysteresis at 0.12~F. These curves further demonstrate that the room temperature cycling amplitude is a decreasing function of gain, until kw > 0.81 where it starts to level off. Figure 7.7 shows the corresponding cycling period behavior. Therefore, it can be stated generally that the room temperature cycling amplitude decreases as the feedback gain of the surface temperature Tw increases and it levels off at kw > 0.8. The cycling period is consistently within practical values i.e. within cycling period values which do not lead to working the mechanical parts of the heating system at a very high frequency. It is to be noted at this point that the practical implementation of a gain different from unity requires more circuitry. Thus, it is recommended for economical reasons to use kw = 1 for the porposed sub-optimum controller.

-147Reviewing the analysis of these graphs indicates a definite improvement in the sub-optimum system over the conventional heating system. This improvement is apparent from the differences in the various room temperature cycling amplitudes. Therefore, it may be concluded that the addition to the conventional heating system of a feedback signal which is proportional to the temperature of the wall surface, improves the performance of that system. If either the sub-optimal or conventional system is to be compared with the optimal system, the basis of comparison must be the defined quadratic performance criterion. It is true by definition that the optimal system developed is the best with respect to this criterion; however, interesting points can be made by analyzing the systems in general. If the conventional heating system is analyzed and compared to the optimum heating system it is found that: (i) the peak to peak variations of the room temperature are much greater for the conventional heating system, when compared to the maximum deviation of the optimal system. (ii) For an outside temperature distrubance the response and adjustment of the optimum heating system is superior to the corresponding response of the conventional heating system. In the optimum system, the temperature begins to fall gradually (due to an outside temperature drop) until it deviates to -O.15~F. Then within about 5 minutes it tends to remain to within 0.1F or less from the otheriginal value. The conventional heating system, on the other hand, begins to oscillate. The rates of increase and decrease in the room temperature depend on the

-148thermostat time constant, and thermostat hysteresis, as explained in Chapter 4. They also depend on the nature of the distrubance. (iii) The rate of change in the room temperature is greater when using the conventional controller thus causing the conventional heating system to be less comfortable. From the analog computer results on the sub-optimal system, it is evident that this system lies in-between the conventional controller and the optimum controller as far as general characteristics are concerned. (i) The room temperature variations are lower than in the conventional system, but higher than in the optimal system; however, the maximum room temperature variations are closer to the optimum system. (ii) The sub-optimum system oscillates between on and off positions when a disturbance is applied, a characteristic of the conventional system. This is anticipated since an on-off furnace system is employed in the sub-optimum system. (iii) The sub-optimum controller tends in general to have a smaller cycling period; however, these values are by no means near to operating the heating system at a high mechanical frequency, which leads to excessive wear of the gas-valve. 7.3 The Practical Implementation of the Sub-Optimal Controller It is evident that the financial cost of the optimal control of a gas-fired forced-air domestic heating system makes impossible the implementation of such a system. The optimal controller devised in this study requires a heat exchanger that is equipped with a continuous modulation system. Such furnace accessories will cost the consumer

-149approximately $200.* A solid state modulation controller is currently produced by the Maxitrol Company. It is used only for several specialized installations requiring rigid performance. The optimal controller also needs four sensors, and the arrangement for feedback gains. Such sensors and gain controls cost at least six times as much as any currently produced thermostat. Therefore, it is evident that the optimal controller, although possible to construct technically, cannot be implemented practically because of current economical considerations. The same cost principle applies to any sub-optimal control system. However, the sub-optimal system suggested in this study was selected on the basis of performance and capital cost. The concepts derived in Chapter 4 indicate that a preferred on-off heating system is one possessing a small thermostat time constant, a small thermostat hysteresis value, and a small heat exchanger time constant. The time constant of the heat exchanger depends primarily on its size, which in turn is determined by the size of the system to be heated. The time constant of the thermostat on the other hand, depends primarily on the properties of the bimetal and the speed of its response, since it is the sensing element. One method of improving the time constant of the sensing device is to employ a thermistor. Thermistors vary in response according to quality which is directly proportional to cost. For example, thermistors of the type Keystone RL21FIT have about a 20 second time constant, and are priced at a very reasonable rate. Such a time constant is much smaller than the time constant of any bimetallic thermostat. The amount of the hysteresis of a bimetallic thermostat depends primarily A figure given by the Maxitrol Company, Southfield, Michigan.

-150on the properties of the relay or the switching device used, Physically, a switching device or relay will produce some hysteresis. However, this hysteresis may be minimized by increasing the sensitivity of the system. The sensitivity of such a thermostat system can be increased by amplifying the signal that is proportional to the room temperature deviation. Therefore, taking these ideas into consideration along with the fact that the addition of a feedback signal that is proportional to Tw improves the performance of the heating system, the sub-optimal controller can be implemented practically at a reasonable cost. The idea of this sub-optimal controller is sketched in Figure 7.10o The temperatures TR and Tw are sensed by thermistors in the sub-optimal controller, rather Power Supply The rm is-te r i ElectricaR + Transr K(TR + Tw) Bridge bTR 6 Transistor Relay ~ T -Bridge "s ~mlfe ~ Relay, Amplifier The rmLster~ l I'I -TRo Gas Valve Reference Reference Figure 7o10. Block Diagram of the Sub-Optimal Controllero

D.C. Power Supply P Tl cm. Rl ^ ~~~~~CI I R2 C A.C. I ^ C R 2 I ICRi Electrical Bridge I20x4 K1~^ - ~~ ~ ~- ~- I I" < t Xl X2 Transistor ^ j~~-I —~ ~ ~ ~ - ~ ~ ~ ~ ~1 ^~ ~R7 i Amplifier i Relaj I Bas-i I I_ Valves I 1XVCircuiti Gas Val ve I i \ R8 ^ ^y \ SCRi I \ ______I Figure 7.11. The Solid State Sub-Optimal Controller.

-152than using a conventional bimetal element. These signals are compared with reference signals TR and Tw, by using an electrical bridge circuit. The difference signal is then amplified by a one stage transistor amplifier. The output of the amplifier actuates a relay. This relay then operates the gas valve of the furnace in an on-off manner. The electrical bridge, the transistor amplifier, and the relay circuits all receive power from the supply through a transformer which must be included. Figure 7.11 shows the detailed solid state sub-optimal controller. The feedback gain, kw, for the surface wall temperature, Tw, is selected as unity in this design. This selection allows for the use of fewer components than if a system possessing non-unity gain was used, and at the same time, it renders considerable improvement to the system. This fact was demonstrated in the previous section. The room temperature reference TR -- is set by the resistor R3. Tw, the.wall surface' temperature reference, is set by R4 Thermistor R5 senses the actual room temperature TR, and thermistor Ro senses the actual surface wall temperature Tw. The output of the bridge circuit which has a signal proportional to (STR + KTw), feeds into a transistor amplifier Q1. The output of this amplifier operates the relay K1, which in sequence turns the gas valve on or off. Preliminary experiments on the sub-optimal controller revealed a 75% reduction in the time constant from the average existing bimetallic thermostat. The hysteresis of this sub-optimal controller, had a value of 0.18~F*. The best bimetallic thermostat available showed approximately 1 F hysteresis in the same test. These values were measured using the standard chamber of the Maxitrol Company, Southfield, Michigan.

It is believed that a sub-optimal controller possessing these characteristics given above including its estimated price may be a significant competitor in the heating industry field. Summary In this chapter a sub-optimal controller for the gas-fired, forced-air heating system has been presented which has improved performance compared with that of a conventional heating system, and which also has an economical practical implementation. This system was simulated on the analog computer and its performance compared to the optimal system and the conventional system. After showing that the sub-optimal controller offers a substantial improvement over existing systems, a practical implementation was presented and shown to overcome the financial problems which beset a truely optimal system.

CHAPTER VIII CONCLUSION The problem of analyzing and optimizing the forced-air gasfired domestic heating system using modern control theory was treated in this dissertation. A mathematical model of the domestic heating system which is acceptable for engineering and for simulation purposes was developed. The model was accurate enough for predicting the effects of parametric changes in system response, and for mathematical analysis in the valid operating region of temperatures. Typical heating system parameter variations were studied using this simulation model. Concepts concerning these system parameters were derived from this analysis and should be of considerable assistance to heating system designerso An optimal heating system for a defined integral quadratic cost function has been developed which incorporates the main objective of minimizing room temperature variations. The optimal control was shown to have the desirable property of providing additional feedback loops to account for disturbances in the system. The feedback portions of the optimum control heating system were also shown to be time-invariant, a characteristic which is advantageous in practice. Parameters of the optimum controller were determined through the use of the Control System Algorithm Program (CSAP) on the 7090 digital computer at the University of Michigan. The resulting optimum heating system was then simulated on the electronic analog computer in order to study its performance in comparison with other systems, namely, (a) the conventional heating system and (b) the sub-optimal system. For the heating system -154

-155under consideration, the optimum controller appeared to drive the room temperature to within 0.1~ F. These results were shown to be superior to the corresponding conventional heating system. The optimum heating system represented an optimum from the theoretical point of view for the configuration and cost function selected. Therefore it represents an upper bound or standard with which conventional.or sub-optimal systems may be compared. However, for some specialized installations possessing rigid performance standards, it may be feasible to utilize a system such as the optimal. Primarily due to economical considerations, it appeared that valuable information might be gained by investigating sub-optimal adaptations of this type of optimum system. In Chapter VII, a sub-optimum control system was proposed that possesses some characteristics of both the optimum and conventional systems.and which is economical. The main difference between the optimal and the sub-optimal controllers is that for the optimal system a modulated heat exchanger or one whose average flame temperature can be varied continuously over wide units is requiredo On the other hand the sub-optimal system used a heat exchanger with onoff characteristics. This sub-optimum system was simulated on the electronic analog computer and compared with the corresponding conventional and optimum systems. The results showed a significant degree of improvement could be obtained using the proposed sub-optimum controller rather than a conventional controller. The practical implementation of the sub-optimum controller was treated in Chapter VII. Since it was shown in Chapter IV that a

-156good on-off heating system would be one that possessed a small thermostat time constant and small hysteresis, a practical controller which possesses small hysteresis and small time constant was devised. Utilizing a controller possessing these characteristics and with ideas devised from the optimal system model a system was practically implemented, and a-solid state controller was constructed. Based upon estimated prices it is apparent that a sub-optimal controller could be made a significant competitor in the heating system industry.

APPENDIX A TIME DELAY CALCULATIONS An element volume of air coming in the room will travel the jet region in T1 minutes, region 2 in t2 minutes and region 3 in 3 minutes. Region \. 2 0 region HRegon 3 co x ~ w. dW Figure A-1. Configuratidn Showing the Different Regions. A-1 Time Delay in Region 1 - X dx 0V.~ 1 o vx But Vx = E vx E x then 1 1.4x v x 3 A ovo T o.8 o H 1.4x x=,15 minutes 0 \/AoV o 1 3 A A-2 Time Delay in Region 2 As the air jet enters the second region it will possess an average velocity vx o It is expected that the air will rise for a short -157

-158time, reflect from the ceiling and then move into region 3. If we assume a linear variation for the velocity in this region then an estimate for the time delay can be obtained. In Chapter II, it was shown that: 3 TA v v = ft/min. X l4x and assuming continuity of flow in the second region then: V2 =. xAx ft/min. wL - AX where v2 average velocity of air entering region 3 ft/min. (1/2) AX = area of jet as air enters region 2. If T2 = time delay in region 2, then:, HX m m min. 2 l/2(vX + v2) A-3 Time Delay in Region 3 A rough estimate of the time delay in this region can be obtained if we assume that the air travels with a uniform velocity v2 which was determined above. If we denote the time delay in this region by T3 then: X = X- m min. 5 v2 Comparing now the three time delays it is found that: T1 - AX T7 2(wL - Ax) 2 A= T3 2wL

-159and since the cross-section area of the jet AX is very small compared to twice the cross sectional area of the room space (about 1% or less), therefore it follows that: TI << T; and T2 3 T

APPENDIX B BASICS OF THERMAL CIRCUITS The idea of a thermal circuit rests on the fundamental similarity between the flow of heat within a rigid body and that of a charge in a non-inductive electric circuit. Conservation of the scalar quantity charge corresponds to conservation of heat. The scalar point function, electrical potential, in electrical field theory corresponds to the scalar point function in heat flow temperature. Ohm's law corresponds to Fourier's law, The concept of electrical capacity of a conductor corresponds to the concept of thermal capacity of a portion of mass. Similarly the concept of electrical conductance of a conductor corresponds to the concept of thermal conductance of a portion of mass. B-1 Flow of Heat in Solids: If the temperature of a body varies throughout its volume, heat will flow from the sections of higher temperature to the sections of lower temperature. Let T(x, y, z, t) be the temperature at a point (x, y, z) within a body at time t. Let A be a plane or curved surface passing through that point and let the coordinate n represent directed distance along the normal to A at the point. Then the flux of heat per unit area per unit time transferred across A at the point by conduction is given by the formula: ~(x,yz,t) = - K T(x,y,z,t) dn Where ~ is the flux of heat which is the rate of flow of heat' and K is the coefficient of thermal conductivity of the material. This formula means that the rate of flow of heat across any face of the volume -160

-161will be proportional to the temperature gradient normal to the face. The negative sign indicates that the heat flow is considered to be positive in the direction in which the temperature decreases. This is the imperical law or the basic postulate for the conduction of heat in solids. If the temperature is only a function of x and t, i.e. T(x,y,z,t) = T(x,t) then the flux in the direction x is given by: (x,t) = -K aT(x,t) ax as 6 aT _ = - Kat ^a-rx where a1 Consider now the element of solid shown in Figure B-1 / (x + Ax) Heat flux Heat f/ / F1moHeat flux by Oen/ Figure B-l. An Element of Solid with Heat Flow in One Direction.

-162The rate of flow of heat through the area by bz at x is given by aQ Ox = at The rate of flow of heat through the area by 5z at x + Ax is given by (X+Ax) = aQ x (Q) bx Therefore the net heat gained by the element is given by: - [Q+ A (Q) bx ] at at ax at or a (I') bx ax at This is equal to the mass of the element multiplied by the rate of rise of its temperature multiplied by the specific heat of the material: aT P 5x by z Cp. a dt Thus, () Cp by 6z T (B-) dx dt at and -K by -z (B-2) at ax where p is the density of the body and Cp is its specific heat. It is to be noted that (dQ)/(-t) in Equation (B-2) is the total heat flux over the area by z. From Equations (B-i) and (B-2) the following equations are obtained: aT- K dT (B-3) at pCp ax2

-163and a2a = K a (a ) (B-4) at2 pcp at -ax It is to be noted that similar expressions could be written for the other coordinates. In a unidirectional flow, the area will not be restriced to 5z by, since it will be constant. Assume it is constant and let it be A. Equations (B-l), and (B-2) then become: - = KA (B-5) at ax - ( Q ) A - Cp (B-6) ax at at B-2 Flow of Charge in a Noninductive Transmission Line: By applying Kirchhoff's laws to a uniform transmission line constructed to have a resistance R and a short capacity C per unit length,it follows that: i = v (B-7) R ax = - (B-8) ax at ai ~+ xx i IR 5x ax _ 6x Figurex B-2. Transmission Line Lumped Model. Figure B-2. Transmission Line Lumped Model,

-164where i is the current, and v is the voltage, as defined in the circuit of Figure B-2. Equations (B-7), and (B-8) are derived as follows: av Voltage change along length 5x = v - [v + x ] ax ax Using Ohm's Law, therefore: R. x i = x ax or I vy 1i = 1av (B-9) R 6x Change in current flow per unit length i - [i + Sx] ax ax 6x And since the current is the rate of flow of charge, therefore: C ax (t) =- a-8x or = -aC (B-bO) ax at From Equations (B-9) and (B-10) the following equations are obtained: av = av (B-ll) at RC ax2 1i = 1 a i (B-12) at RC ax2 Note that i = (dq)/(dt), where q is the quantity of electricity.

-165B-3 The Analogy between Thermal and Electrical Constants and Variables: Since Equations (B-3), (B-4), (B-5), and (B-6) are respectively of the same form as Equations (B-li), (B-12), (B-7), and (B-8), the following analogy can be maintained and thus an electrical respresentation may be made of any heat problem. Heat Problem Electrical Problem Thermal Conductivity x Area Conductance/unit length of trans(KA) mission line 1/R Area x density x Specific heat Capacitance per unit length of i.e. thermal capacity per unit transmission line length (ApCp) (C) Quantity of Heat Quantity of Electricity (Q) (q) Flow of Heat Current (i) Temperature Difference (T) Potential difference across transmission line (v) B-4 Scale Problems The values of the electrical quantities representing either variables (T, ('Q)/(6t), x and t) or constants (K, P and Cp) in the thermal problem may be scaled by any factor. However, changes have to be done in order to preserve the equalities of Equations (B-7) and (B-8). If the value of the capacity is multiplied by a factor n both the equalities (B-7) and (B-8) are preserved merely by modifying the time scale. Equation(B-7) remains unchanged and Equation (B-8) becomes: x = -nC ~x nan~t

-166If the value of the resistance is changed, the equalities cannot be preserved by altering only one parameter. One way is to modify i and t. If the scale of i is divided by n, Equations (B-7) and (B-8) become i 1 av m mR ax 1 ai av m ax mat It is noted that since R and C occur in Equations (B-li) and (B-12) only as the product RC, the result will not be affected by changing R and C individually, provided their product is not changed. Also since Equation (B-ll)is homogenous in v, v may be replaced by 1v without affecting anything else. By suitable choice of scale factors m, n, & the electric model may be operated at convenient voltages and transit time intervals and built with feasible magnitudes of resistance and capacitance. A heat process actually taking hours or days may be condensed to a few minutes in the experiments. On the other hand, a short process requiring only fractions of a second may be stretched in the model so as to last several minutes.

APPENDIX C HEAT CAPACITY ESTIMATION OF THE ROOM BOUNDARY In this appendix an experimental check is made for the estimation of the room boundary equivalent capacity. The room boundary parameters are calculated first. For every layer of the boundary, the heat capacity and thermal resistance are calculated. Then these capaciters and resistors are connected to form an RC network representing the boundary. For region two, which faces the lower part of the wall and the floor, the boundary model is two RC networks in parallel. The same applies for region three. The electrical circuits representing the models O'F C~ 1 T hRC Netgeerto ic ee csh nwork of 0 1) C) of the output is observed by an oscilloscope to measure the time constant of the model. Knowing now this time constant, and the thermal resistance of the one T-network model calculated in Appendix D, the equivalent capacity may be estimated. -167

-168C-1 Estimation of the Capacity of Region Two-Boundary: Region two faces: an area of 91 ft2 of the wall and the floor. Figure C-2 illustrates the construction and electrical model of the wall. Using Table C-I, and the formulas: Heat Capacity C = density px specific heat Cpx volume V; and TABLE C-I PROPERTIES OF THE BUILDING MATERIALS Material p lb/c.f. Cp Btu/lb. ~F k Btu hr ft? ~F plywood.84.57.085 Aluminum 2.8.266.122 Air.0808.2375.016 Fiber glass.022 Tempered Hardboard.71.65.125 Wood.65.57 0.085 thickness x Thermal Resistance R = -thickness x, which were Thermal Conductivity k x Area A defined in Appenidx B, the thermal resistances and capacitances are calculated for the section facing region 2'. R1 = 5 ~F/(Btu/sec) 1 = 1.36 Btu/~F R2 = 69 ~F/(Btu/sec) C2 = 0 Btu/~F R3 =. O F/(Btu/sec) C3 = 0.12 Btu/~F R4 = 118 ~F/(Btu/sec) C4 = 0.18 Btu/~F R =- 0 ~F/(Btu/sec) C5 = 0.15 Btu/~F R6 = 1.67 ~F/(Btu/sec) C6 =. 45 Btu/~F

-169Similarly for the floor: R1 = 32.4 ~F/(Btu/sec) C1 = 15.4 Btu/~F R2 = 255 ~F/(Btu/sec) C2 = 0.884 Btu/OF R3 = 13.3 F/(Btu/sec) C3 = 6.38 Btu/~F These two networks are connected in parallel to form the electrical model of the boundary facing region two. A measurement for the time constant of this boundary was made using a square wave generator and an oscilloscope. A time scaling was introduced in the electrical analog so that 1 sec on the electrical model corresponds to 106 seconds in real time (Refer to Appendix B for scale problems). The measured time constant was 2.1 milliseconds which corresponds to 2100 seconds in real time. If the response of such a model is to be approximated with a one RC T-network, whose resistance is the thermal resistance of the region two boundary, then the value of the capacity of the boundary would be estimated as: Time Constant measured Thermal resistance of Boundary + Impedance of generator 2100 - 12.45 Btu/OF 1.975 x 60 + 50 C-2 Estimation of the Heat Capacity of Region Three-Boundary The same method is employed to estimate the heat capacity of the boundary facing region three. This region faces a 354.2 ft2 of the wall shown in Figure C-2, and the ceiling shown in Figure C-4. The RC parameters of the wall portion and the ceiling are:

-170 - I I I I. I 3/8" 1" 0.025" I 5/ 4" 0.032" 1/8 5/8' Aluminum I Tempered Hardb Ixd Plywood I Fiberglass Aluminm I Air Space Aluminum Tempered Hardba R1 B1 1R2 R2 R3 1 R4 R4 R5 R5 I R6 R C y2 y_ i TC I T - T,, Figure C-2. Construction and Electrical Model of the Wall. 3-5/8" I 4 1" Wood Air \ Plywood R R R R R R | 1 1 | 2 R2 I 3 32 3 T 1l I I! Figure C-3. Construction and Electrical Model of the Floor. 1 3/8" 1 1".032 " 3" I 0.032" 1 1/8" I Plywood Fiberglass \ Aluminum I Fiberglass I Aluminum Tempered i ~~~~~I I ~Hardboard Ri R1 R2 R2 R3 R3 R4 R4 I R5 R5 R6 R6 I I FI I Igr i Figure C-4. Construction and Electrical Model for the Ceiling.

-171For the part of wall facing region 3: 1 = 1.28 ~F/(Btu/sec) C1 = 5.31 Btu/OF R2 = 17.7 ~F/(Btu/sec) C2 = 0 Btu/~F R3 = 0 OF/(Btu/sec) 3 = 0.468 Btu/~F R4 = 0.303 ~F/(Btu/sec) C4 = 0.702 Btu/~F R5 = 0 F/(Btu/sec) 5 = 0.585 Btu/OF R6 = 0.428 ~F/(Btu/sec) C6 = 1.75 Btu/OF For the ceiling: R1 = 1.12 OF/(Btu/sec) C1 = 0.69 Btu/~F R2 0 ~F/(Btu/sec) C2 = 0.233 Btu/~F R3 = 140 F/(Btu/sec) C5 - 0 Btu/OF R4 0 ~F/(Btu/sec) C4 = 0.233 Btu/OF R5 46 ~F/(Btu/sec) C5 - 0 Btu/OF R6 = 3.4 OF/(Btu/sec) C6 = 2.06 Btu/~F Time constant measured -:= 2.2 millisecond Estimated Capacity 2.2 x 10 0.66 x 60 + 50 =24.5 Btu/OF The capacities determined in this appendix using the electrical analog are within approximately 7% of the corresponding values measured experimentally in Chapter III. It is to be noted that the accuracy of the electrical model depends on the accuracy of the data of parameters for building materials. There is some uncertainty in these data because

-172of the porosity of the materials and uncertainty of their water content. Thus, it follows that there is also some uncertainty in the estimated values of capacities using the electrical analog method. On the other hand, an error results when using the experimental measurement (used in Chapter III) due to the accuracy of the measuring equipment. Therefore, it appears that there is some uncertainty in both methods. The values estimated experimentally were used for the mathematical model in Chapter II, and it resulted in an average percentage error of less than five when analytical and experimental results were compared.

APPENDIX D EXPERIMENTAL MODEL CALCULATIONS AND TEMPERATURE PLOTS D-1 Volume Calculations From Figure D-l the volume V, of the jet region is calculated: V1 = (1.94)2(7.75 + tan 11~) = 30.7 ft3 11~ 7\ 75 - 3.88- Figure D-l. Dimensions of Region 1. Volume of region 2 = V2 = 11.76 x 11.76 x 1.94 = 270 ft3 Volume of region 3 = V3 = 11.76 x 11.76 x 7.75 - 270 = 10493 ft3 D-2 Thermal Conductivity Calculations: In order to calculate the equivalent thermal conductivity of the walls, the ceiling, and the floor of the experimental room model the following equation for a multi layer material is used:(3) x x Xi (D-1) i=l ki -173

-174where, Xi = thickness of layer i, ki = thermal conductivity of layer i, x = total thickness under consideration k = equivalent thermal conductivity of the multi-layer material. Applying now Formula (D-l) to the wall, which is composed of six layers, then: Xw 6 xi kw i=l ki where, 6 x, = total thickness of wall = Z xi = 2.807", k = equivalent thermal conductivity of the wall. w Using Table D-I, kw is calculated: kw = 0.0238 Btu/ft ~F hr. TABLE D-I a. WALL PROPERTIES Thickness Thermal conductivity Material inches Btu/ft ~F hr Plywood 3/8 0.125 Fiber Glass 1 0.024 Aluminum 0.025 118 Air 1.25 0.0175 Aluminum 0.032 118 Tempered Hardboard 1/8 0.125

-175TABLE D-I (CONT D) b. CEILING PROPERTIES Thickness Thermal Conductivity Material inches Btu/ft ~F hr Plywood 3/8 0.125 Fiber Glass 1 0.024 Aluminum 0.032 118 Fiber Glass 4 0.024 Aluminum 0.032 118 Tempered Hardboard 1/8 0.125 c. FLOOR PROPERTIES Thickness Thermal conductivity Material inches Btu/ft ~F hr Plywood 1.5 0.125 Air 4 0.0175 Wood 3-5/8 0.125 For the ceiling: kc 6 Xi kc i=l ki where, Xc = total thickness of ceiling = 4.564", kc = equivalent thermal conductivity of ceiling.

-176Therefore by using Table D-I: kc = 0.0267 Btu/ft OF hr For the floor: _ = z ^ kf i=l ki where, Xf = total thickness of floor = 9.125", kf = equivalent thermal conductivity of floor Hence, using Table D-I: kf = 0.0338 Btu/ft OF hr D-3 Thermal Resistance Calculations In this section two main parameters for the experimental model are calculated: R2 and R3 the thermal resistances of the boundaries surrounding regions 2 and 3 respectively. R2 is the parallel combination of two resistances: RB which is the thermal resistance.2 of the part of the wall facing region 2 and Rf which is the thermal resistance of the floor. Similarly R3 is the parallel combination of Rw the thermal resistance of the part of the wall facing region 3 and Rc the thermal resistance of the ceiling. By using the analogy between heat and electrical problems explained in Appendix B, then:. xf "v Rf = -~ ) oRw2 = kfAf kwAw

-177and xc xw Rc= kcAc' Rw3 krA where AfAw2,A and Aw3 are the areas of the floor, the part of wall facing region 2, the ceiling, and the part of wall facing region 3 respectively. Af 134 ft2, A2 = 91 ft2 Ac = 134 ft2, Aw3 = 354.2 ft2 Calculating, then: Rf = 10.1 OF/(Btu/min) Rw2 = 6.5 OF/(Btu/min) Rc = 6.38 OF/(Btu/min) R3 = 1.66 OF/(Btu/min) Also =1 =f PW2 = 1,975 ~F/(Btu/min) 2 Rf+Rw2 and R Rw R 1 RcR o0.66 ~F/(Btu/min) 2 Rc+Rw

-F ~A 100 - 90 so - 80 70 - 60 50 - 401 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 TIME, min Figure D-2. Temperatures of Points in Region 1.

hO 90 F 80 so 70 60 500 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 TIME, min Figure D-3. Temperatures of Points in Region 2.

OF 100 G 90 - 80 - 70 60 4 1 I 30 40 50 I0 1 140 1 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 TIME, min Figare D-4 Temperatures of Points in Region 3

*F 10o 90 80 K L so ~~ -^ ^ K 70 60 50 401- 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 TIME, min Figure D-5. Surface Temperatures of Points Facing Region 2.

F 100 90 80 M_ N ~~~~~60 ^^^^~~~~~~0 70 60 50 40I 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 TIME, min Figure D-6. Surface Temperatures of Points Facing Region 5.

REFERENCES 1. "Distribution of Air Within a Room for Year-Round Air Conditioning," Parts I and II, Harold E. Straub, Stanly F. Gilman, Seichi Konzo, and Michael Ming Chen. University of Illinois Engineering Experiment Station Bulletin Nos. 435 and 442. 2. American Society of Heating and Refrigerating and Air-Conditioning Engineers, Guide and data book, 1963, Chapter 11. 3. "The Study of the thermal behavior of Structures by Electric Analogy, G. Brunard, British Journal of Applied Physics, 1950. 4. "Room Temperature Response to a sudden heat disturbance input," by R.O. Zermuehlen, H.L. Harrison, ASHRAE 1965. 5. "Dynamic characteristics of water to air cross flow heat exchangers," by J.R. Gartner, H.L. Harrison, ASHRAE 1965. 6. Warm Air Heating and Air-Conditioning Digest," Indoor Comfort." 7. American Society of Heating, Refrigerating and Air Conditioning Engineers, Guide and data book, p. 152. 8. "Time modulations control increases comfort," Part I and II, by Ralf Allen, American Artisan, July 1950. 9. "Cooling Load from Pretabulated Impedances," by Harry Buchberg, Heating, Piping and Air Conditioning, February 1958. 10. "Cooling Load from Thermal Network Solutions," by Harry Buchberg, Heating, Piping and Air Conditioning, October, 1957. 11. "Transfer Admittance Functions," by W.B. Drake, ASHRAE, 1959. 12. "The Optimal Control of Gas-Fired Warm Air Heating Systems," by L.F. Kazda, G. Casserly, A. Eltimsahy and R. Spooner. AGA Final Report Contract DO-14-GU, University of Michigan. 13. Heat Transmission, by W.H. McAdams, McGraw-Hill, 1942. 14. Optimization Theory and the Design of Feedback Control Systems, C.W. Merriam III, McGraw-Hill, 1964. 15. "Air Velocities in two parallel Ventilating jets," by A. Koestel and J.B. Austin, Jr., Transactions of American Society of Heating and Air-Conditioning Engineers,Vol. 62, 1956. 16. Introduction to heat transfer, by Brown and Marco, McGraw-Hill, 1958. -183 -

-18417. "The Method of Analogue Computing used at the Heating and Ventilating Research Association Laboratories," J.I.H.V.E., July, 1964. 18. Introduction to Electronic analog computers, by J. Warfield, Prentice Hall, Inc., 1959. 19. "Two Position discontinuous temperature control in Electric Space Heating and Cooling processes," W.K. Roots, J.M. Nightingale, IEEE Transactions on Application and Industry, Jan., 1964, pp. 17-27, Part I. 20. "A Survey of discontinuous temperature control method for electric space heating and cooling processes, " W.K. Roots, J.M. Nightingale, IEEE Transactions on Application and Industry, Jan., 1964, pp. 1-12. 21. "Two Position discontinuous temperature control in Electric Space Heating and Cooling Processes," W.K. Roots, J.M. Nightingale, IEEE Transactions on Application and Industry, Jan., 1964, pp. 27-38, Part II. 22. "Multiposition discontinuous Temperature Control in Electric Space Heating Process," W.K. Roots, J.M. Nightingale, IEEE Transactions on Aplications and Industry, Jan., 1964. 23. Modern Control Theory, J. Tou, McGraw-Hill, 1964. 24. Linear Systems, R. Shwarz and B. Friedland, McGraw-Hill, 1965. 25. Analog and Digital Computer Technology, N. Scott, McGraw-Hill, 1960. 26. "Comparative Performance of Year-Around Systems Used in Air-Conditioning Research Residence No. 2," by J.K. Wright, D.R. Bahnfleth, E.J. Brown, University of Illinois, Bulletin 465. 27. Network Analysis and Synthesis, by F.F. Kuo, Wiley, 1962. 28. "Supply Outlet Locations for Basement Heating," by J.R. Wright, and D.R. Bahnfleth, ASHAE Journal section, Heating, Piping and Air-Conditioning, August, 1958. 29. "Investigation of Summer Cooling in the Warm-Air Heating Research Residence, " University of Illinois Bulletin No. 36. 30. Heat Transfer for Electrical Engineers, by A.D. Moore, Wahr Pub. Co., Ann Arbor, Michigan, 1949. 31. "Heat Emission Characteristics of Warm-Air Perimeter Heating Ducts," University of Illinois, Bulletin no. 412. 32.'"Dynamic Programming and the Calculus of Variations," Stuart E. Dreyfus, Journal of Mathematical Analysis and Applications 1, pp. 228-239., 1960.

-18533. "Optimization Theory and the design of Feedback Control Systems, C.W. Merriam III, McGraw-Hill, 1964, Chapter 4. 34. "Fundamentals of Heat Transfer in Domestic Gas Furnaces," Research Bulletin 63, American Gas Association Laboratories. 35. "Control Systems Analysis Program II, "(CSAP II), User's Manuel. 36. American Society of Heating, Refrigerating and Air-Conditioning Engineers Guide and Data Book, 1963, p. 185.

UNIVERSITY OF MICHIGAN 3 9015 02826 6735