THE U N I V E R S I T Y OF M I C H I GAN COLLEGE OF ENGINEERING Department of Mechanical Engineering Student Project Reports INVESTIGATION OF DESIGN MEANS FOR HOME LAUNDRY APPLIANCES Jack P. Smith Kenneth J. Timmer ORA Project 371550 supported by: WHIRLPOOL CORPORATION BENTON HARBOR, MICHIGAN administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR July 1971

TABLE OF CONTENTS Page LIST OF FIGURES v Simulation of the Balancing Mechanism on the Combination Washer-Dryer by Jack P. Smith 1 1. GOALS OF THE SIMULATION PROJECT 3 2. DEVELOPMENT OF A MODEL 3 A. Description of First Model 3 B. Equations of Motion 6 C. Revised Model 10 3. SOLUTION TECHNIQUE 15 A. State Space 13 B. Analytical Method 19 4. COMPUTER SIMULATION 23 A. Subroutine CALL 28 B. Subroutine TANK 28 C. Subroutine FIND 28 D. Subroutine OUTP1 31 E. Subroutine OUTP2 31 F. Input Variables 33 G. Output Control 33 H. Output Columns 34 I. Output 35 J. How to Use the Program 45 K. Results and Conclusions 45 REFERENCES 48 iii

TABLE OF CONTENTS (Concluded) Page A Computer Simulation of a Ventless Clothes Dryer Using a Refrigeration Cycle by Kenneth J. Timmer 49 1. OBJECTIVES 51 2. RESULTS 51 3. BACKGROUND 51 4. BASIC EQUATIONS 54 5. ANALYSIS OF THE REFRIGERATION CYCLE 56 6. DISCUSSION OF MAIN PROGRAM 58 7. DISCUSSION OF PSYCKT SUBROUTINE 66 8. DISCUSSION OF GRAPH SUBROUTINE 73 9. DISCUSSION OF DRYER SUBROUTINE 75 10. EXPERIMENTAL RESULTS 77 APPENDIX 85 iv

LIST OF FIGURES Figure Page 1. Pictorial diagram of machine. 4 2. First model. 5 5. Three masses lumped into one effective mass and eccentricity. 9 4. Final model. 11 5. Ballast water flow rate vs. tank displacement. 14 6. Flow diagram for state space solution. 15 7. Manner in which mass accumulates in ballast tank #2. 16 8. Analytical approach towards the solution. 20 9. Effect of damping on time response. 22 10. Effect of damping on phase angle and magnification factor. 24 11. Main program flow diagram. 25 12. Subroutine flow diagrams. 29 13. Subroutine flow diagram. 30 14. Subroutine flow diagram. 32 15. Time to attain balanced spin-data from program output. 46 16. Time to attain balanced spin-at a different flow rate. Data from program output. 47 1. Schematic diagram of ventless dryer system. 52 2. Comparison chart of power and time required to dry a normal clothes load (9 lb dry weight, 85 initial moisture content). 53 3. Temperature-entropy diagram of refrigeration system. 57 v

LIST OF FIGURES (Concluded) Figure Page 4. Schematic diagram of refrigeration system. 57 5. Clothes load and moisture content vs. drying time. 78 6. Air flow vs. drying time. 80 7. Enthalpy vs. drying time. 81 8. Condenser capacity vs. time. 82 A-1. Schematic diagram of open-loop system. 86 A-2. Compressor curves. 89 A-3. Compressor curves. 90 A-4. Compressor curves. 91 vi

SIMULATION OF THE BALANCING MECHANISM ON THE COMBINATION WASHER-DRYER Jack P. Smith 1

1. GOALS OF THE SIMULATION PROJECT The first step in working out a computer simulation of the self-balancing washer was to formulate a model of the system. It was recognized that this model would involve those aspects of the machine pertinent to the high-speed spin condition when the balance control is in operation. From this an analytical model was to be derived including its equations of motion. Second, it was desired to develop a computer simulation from this model which would, for a given set of input conditions, give the time response and other important data on the system operation. Third, the computer program was to be used as a design tool. During this phase of the project the computer model would be run under different conditions and an attempt made to glean from the results possible design recommendations. Further it was desired to show by example how the program could be used as a design tool. 2. DEVELOPMENT OF A MODEL A. DESCRIPTION OF FIRST MODEL Figure 1 is a diagram representing the combination washer-dryer. As the figure indicates, the tank is supported from the bottom by the base plate attached to a channel frame. The position of the balance control device is also shown. The diagram is the mirror image of the actual configuration. This is so that the cylinder rotates in the counterclockwise direction, thus simplifying the trigonometry. Figure 2 is a diagram of the first model. The suspension of the tank has been represented by the torsional spring on the bottom. The geometric center of the tank is allowed to move in a plane perpendicular to the tank's axis. Either rectilinear or cylindrical coordinates could be used. Here the cylindrical coordinates r and * locate the tank's moving geometric center relative to its static center. The belt which drives the tank is elastic. This effect causes Q= 0 e+ 7 3

Figure 1. Pictorial diagram of machine. 4

Tank 1 \ -— To^<^Tank Rotating k,,/r ) S^ \f^ ^^V^ Cylinder Tank 2 Tank 3 a = r + u S= +r FIRST MODEL' = -r2 Figure 2. First model. 5

which says that Q, the angular velocity of the tank, is equal to e, the steady state velocity, plus A, a superimposed oscillatory velocity resulting from the belt's elasticity. Likewise for the drive take off, a = a + 4 where a is the velocity of the drive pulley. a is the steady state velocity, and p is the oscillatory velocity. The two steady state coordinates will be related by the equation r a= - - r 2 where rl/r2 is the ratio of the two pulleys' radii. Thus this is a 5 degree of freedom system.., y, p, r, and 4 are independent coordinates. a is related to 0 by the above constraint. B. EQUATIONS OF MOTION It would be at best difficult to derive the equations of motion for a system such as this by a Newtonian or force-balance approach. Also a hazard of such an approach is that unexpected cross product terms are easily overlooked. A better method is to use the Lagrange energy formulation. The Lagrange equation is: d (T-V) -(T-V) In this equation T is an expression for the kinetic energy of the system V is the potential energy of the system q. corresponds to a generalized coordinate Q. corresponds to the associated nonconservative force, 1 i = 1,2,3,...,n. First, T and V are expressed in terms of the five coordinates, and then the 6

equation is applied once for each coordinate where the coordinate replaces qi lT = /2[2 m2 + mMp2 + m424 + 2 + +2 + i( +)2+ )2] 11 m2P2 m pP3 4P4 T Pi = r + r. 1 1 -- = + r(J X -*) + (6 =~( Pi = r r + r( ) + (+7)(kxri) ()2 = r + (tr)2 + (6 +)2 r - 2(G+7)r ri sin (i i i + 2V( +y)r ri cos ~. 1 1 where.i is the angle between r and ri and ~i = 8 + 7 + ai - f $ A ( r = r r+ k(kxr) -2.2 2 r2 = r + (r)2 So an expression for T can be written: T = 1/2, mi(r2 + (fr) + ( +) r2 - 2(0+y)r ri sin ( i=l + 2(Y+ y)r ri cos ~i) + M(r2 + ( 2r)2 + I( +)2 + l * + 2 The potential energy of the system, V, is determined by superposition. To start the coordinates are all set to zero then the system is changed such that one coordinate at a time has some positive value. This can be seen vectorially by q. 1 v. = f (Z F.) dq. o1 n The total V then is Z Vi. The resulting equation for V is i=l 7

2 2 k r + k r V = 1/2 k r + r M g sin V + 2 sin2 x T 2 4 + E m.gr(sin(O +y +.) - sin a + V. _ i i belt where Vbelt is the energy stored in the belt. The nonlinear terms in the expressions for kinetic and potential energy show up in the equations of motion. It is very desirable to get equations which are linear since the system is basically linear. If this can be done the equations come within the scope of easier solution techniques. The replacement of r and i, the cylindrical coordinates, by x and y rectagonal coordinates reduces some of the nonlinearity but nonlinear terms still prevail in 0 and 7. By lumping the masses into an effective mass and eccentricity, it was hoped to simplify some of the complexity associated with terms involving the four separate masses. As shown in Figure 3, this effective mass would be 4 M = Z m. e l i=l Its eccentricity would be determined by the expression 4 M-e = Z m. r. i=l This equation would indicate that the M and e are the equivalent in a static situation. But since they are to be used in a dynamic expression, they should also be equivalent with respect to momentum. Thus d - d i - (M.e) = - Zmr dt e) dt m r ee 4,. M e + M e = (mf r + m r.) i=l 1 When it was attempted to prove the above relationships, it was found that there was no way to satisfy both constraint equations simultaneously. 8

MI r4 r? 2 rs Figure 3. Three masses lumped into one effective mass and eccentricity. 9

C. REVISED MODEL At this point some time was spent running the machine and determining its operating characteristics. The drive belt is so stiff that it was decided any superimposed oscillation would be small. The tank was found to oscillate about a point 20 in. below the tank center. Due to the stiffness of the system, this oscillation has an amplitude of only.016 in. when the clutch bleed valve is first touched. This is such a small arc that the motion is essentially horizontal and linear. Vertical motion is virtually nil due to the extreme stiffness in this direction. Thus it was possible to revise the model as shown in Figure 4. This is now a 2 degree of freedom system with x and 0 as coordinates. The equation for this system can be found, neglecting the inertia in the transmission: T = 1/2L m()2 + 2 i=l Pi = (x + r cos(+ +a.))l + (r sin(e +a.)) or Pi = xl + r. i i 1 1 (P.i) = + x 2 r. - 29 r.x sin( +a.) so 4 T = 1/2 - 2 26 rx sin( +a.) + Z m. 2 r + M x + I MTOT wi=lr i= where 4 M M + m. TOT i 1 1 10

r2 44 / ^ K11 \\\ J\Z< r. 1~~~~~1 I~~~~' I —, 4. F —---- ---' i Figure 4. Final model.

The potential energy is V = V + VE x 0 V /2 kx2 x 4 e+a. V0 = Slrmg cos q (k) * dq.(k) u-i i i=l a. 1 4 V = rm.g(sin(0 +.) - sin C.) i=l 4 V = 1/2 kx + Z rm.g(sin(9 + a) - sin a.) 1 1 1 i=l For the x coordinate the equation of motion derived from the Lagrange formula using the above expressions for T and V is 4 4 Or m sin( +.) + 6r - m sin(e+ i i=l i=1 4 4 + r -l m cos(0 +a) + I mx + M x + kx = 0 i. 1 i i TOT i=l i=l The second and fourth terms in this equation represent momentum transfer to the water being added to the ballast tanks. But since the water enters the ballast troughs tangentially and with a velocity not greatly different, the force on the tank due to this mass addition will be small and these terms can be neglected. Dropping these terms and rearranging 4 4 MTT kx= r m cos(9 +a.) + Or mi sin(6+ +a) i=l i=l where the first summation represents the four centrifugal forcing functions and the second summation is the mass angular acceleration forces. This equation is a, second-order differential equation, linear in x, with time varying forcing functions. The mi are of course the three ballast tank masses and the off-balance load mass. The tank masses are going to be changing but in sequence. At any time, one tank will be changing, and the other two will have constant mass. MTOT changes with time also, but not very substantially. 12

3. SOLUTION TECHNIQUE An effort was made to apply the solution techniques available in control theory to the problem. CSAP (Control Systems Analysis Program) was considered but found to be too narrow to apply to this problem. The fact that the ballast tanks change in spurts or pulses necessitates the need for switches or pulsers in the control diagram to let the signal through the proper feed back loop at the proper time. Both Z and P transforms were investigated but lack of information on these techniques hindered a successful development. A. STATE SPACE State space was next investigated and this technique appeared to be applicable. In state space an nth order problem can be broken down into n first order systems. These are arranged in a matrix form and operated on by various matrix operations to yield a solution equation. This equation looks like a matrix form of the Duhamel superposition integral used in vibration study. Figure 5 shows the relationship between x, the tank displacement, and the flow rate into the ballast tanks. This characteristic was found experimentally and it assumes a linear increase in flow as the displacement ineas. In actuality the flow is probably increasing most rapidly when the flow of circular cross section is half uncovered than when it is only slightly uncovered. But in the interest of being able to find a solution this linearity was assumed. Figure 6 is a diagram representing the state space system. In this system the feed forward transfer function is forced by the four centrifugal forces. The unbalance load force is of constant amplitude. The forces from the ballast tanks are of increasing amplitude where the amplitude of any one varies as the integral of x and increases only when its feed-back switch is closed. Figure 7 shows how mass accumulates in the tanks as a monotonically increasing function with time. The state space problem is set up as follows. Let x = x 1 2 - — i + (^1u 13

bJ n 1.0 -J _ -.0695 --------- z.5./.0158 0.01 0.02 0.03 DISPLACEMENT (INCHES) Figure 5. Ballast water flow rate vs. tank displacement. 14

STATE SPACE cos (wt-y4) F mrw u 1T... X s + K X /MT cos(wt-y) A Y, z 2 w r de lay cos 2 t- 2),'pulse width Y2 z2 a 2 period S T, h, h cos(wt- 3) Y3 z3 2r./ Ct, r S T, h, 2h A = - h = T/3 m Figure 6. Flow diagram for state space solution. 15

| Time Mass (Tank I) Time Mass (Tank 2) X- lz=- e- -- Time NO MASS ACCUMULATES IN TANK 3 Figure 7. Manner in which mass accumulates in ballast tank #2. 16

where MT is the total mass of the system. The above equation can be rewritten in matrix form as 0 1 O (x) = (x) + u - 0 %L T M where (x} =1 (X f 2J and 2 u = m4 r cos(ot t-y +y4 + Y y +Y3 cos (ct -y )Z Y2 = cos (ct- 2)z2 y = cos(Wt-7 )z r Ax kT < t < kT + T/3 k = 0,1,2,... ) 0 ) Votherwise r Ax (kT + T/3 < t < kT + 2T/ k = 0,,.... z ( j when 2 0 J otherwise 2) r2 AxA (kT + 2T/3 < t < (k+l)T k = 0,1,2,.... Z7 = ( /when < t 0 J~ w otherwise where A is the linear relationship between displacement and flow rate X For the first time period kT < t < kT + T/3, the equations for x and Z can be combined. 17

0 1 0 0 0 0 |-k cCo(ot-yl) cos(wt-y2) cos (t-y3 ) m4rw cos (ot-74) - 0 -' - --- MT MT MT M M (r) = Aw r 0 0 0 0 (r) + 0 o 0 0 0 0 0 o O 0 0 0 0 where 1 x2 (r) = I (zl 2 The equation for the second time period, when the second ballast tank is being filled or when kT + T/3 < t < kT + 2T/3, will be the same except the Aw2r will occupy the 4,1 position. In the equation for the third period when the third tank is filled, the Aw2r term will occupy the 5,1 position. Letting [B] represent the square matrix and (C) the force vector (r) = [B](r) + (C) The transition matrix [O(t,t )] required in the solution equation is now found t t 1 [$(t t)] = [I] + S [B(T)]dT + J [B(T)] [B(T2)]dT dT1 +... to t t o Because the problem has time varying coefficients this series form for the transition matrix I is necessary. In a problem with constant coefficients it is possible to obtain the transition matrix in closed form. Once the transition matrix is determined the solution equation is 18

t (r(t)} = [$(t,t )][r(t))} + f [9(t,T)]{C(t)}dT t 0 This equation looks a lot like the Duhamel superposition integral used in vibration study. It contains a transient term dependent upon the initial conditions and a steady state portion, the second term. The clear advantage of using state space is that once the transition matrix [0] and J[O](C)dT are determined, the response of the system is easily obtained at any time t given the initial conditions at time to. In this particular problem there will be three different solution equations, one of which will correspond to each tank. And that equation will be valid while its corresponding tank is adjacent to the flow nozzle. So a time increment of 1/3 of a revolution or 120~ is possible. In attempting to obtain the solution equation, the first five terms of the transition matrix were determined. As each successive term is the integral of the product of the [B] matrix times the preceding term, and these matrices are 5 by 5, the elements of the fourth and fifth terms' matrix become very involved expressions. It was hoped that five terms would be sufficient to assure convergence, but due to the complexity of the terms it was very difficult to evaluate by hand. As a result, a computer program was developed based on this solution method. The results obtained after the program was running indicated a basic error in formulation. It turned out that in the second term of the solution equation the integral was taken over t rather than to. This change makes a big difference in the steady state term. Also there was question as to whether five terms assured convergence of the transition matrix. B. ANALYTICAL METHOD At this point it was decided to take smaller time increments, assume coefficients are constant over that increment and solve the problem analytically. This scheme is shown in Figure 8. As indicated in Figure 5, there is a maximum flow rate point past which an increase in X does not increase m. Also for negative X the flow is zero. Thus in Figure 8 there are two solution schemes. When X is positive, small time increments are taken and the mass change is determined by integrating the flow rate over the time increment. This change is added to the appropriate ballast tank, the initial conditions are reset and another time increment is taken. As the displacement may increase past the point at which mmax has been attained, it is necessary to determine the time at which this occurs, marked T1 on Figure 8. Likewise T2 must be determined. When X becomes negative there is no flow rate so no mass changes, and the time increments can be increased. In Figure 8, 30~ time increments are shown. It is necessary to determine the time at which X = 0 so the appropriate solution technique can be applied. 19

thf x1 MAX Discre Zero-finding operations Analytical Solution T -., -u —-- — Time TI T2 Analytical Solution X=eE'~Unt (Aicos Wdt+A2Sin Wd t)++ K Sin(wt- i - c -) Figure 8. Analytical approach towards the solution.

Originally the system was considered to be undamped. When the program based on the undamped solution was run, a time response as shown in Figure 9(a) was obtained. The high frequency superimposed signal is due to the transient portion of the solution and in an undamped system this does not attenuate. The real system has a certain amount of damping and does not behave this way. Although the damping in the real system is structural or hysteretic it was found that the addition of a small amount of viscous damping into the model attenuated the transient, "filtering" it out of the steady state signal. This response is shown in Figure 9(b). The equation of motion with damping added becomes 4.2 M X + cx + kx = 0 r Z m. cos(0 +C.) + r Z m. sin(+i+.) TOT i=l i=l The analytical solution to this equation is 2 4 m. rc X = e-nt (A cos wdt+A2 sin wdt) + E K sin(ut +- +.- ) 1 d 2 d. 2i=l 4 m.ro + ZE 1 K sin(t +aY. - i=l k 1 where o = 0 the angular speed and XO = 0 the angular acceleration. From Figure 4 it can be seen that: 1' the position of ballast tank 1 = -240~ a2, the position of ballast tank 2 = 0~ a, the position of ballast tank 3 = -120~ a4, the position of the off balance load is equal to -y4 in other notation, and is variable. Al and A2 are constants of integration which depend upon the initial conditions. 2 i * 4 m.r! 4 m.ro A = X() - Z 1 K sin(-+a. - K sin(a.- 1 k 2 i k 1 i=l k i=l 21

50 RPM 50 RPM 5 lb. unbalance =. 1 1.0 UNDAMPED DAMPED RESPONSE RESPONSE.8.6.4 V).2 "r 1i1 Vi' I.o 01. 2.5 -5 1.0.2.4.6.8 -1. 0 (a) (b) Figure 9. Effect of damping on time response.

_1 r, 4 m. ro A = + X(0) + A n Z - 1 K cos(-+ai- 4) 2 wd n i=l k 4 m. rc ] - Z 1 K cos(a. - 1) k 1 i=l 4 is the phase angle between force and displacement which is present in a damped system. 4) = tan- 2-r l-r where r here is the frequency ratio w/cn, not to be confused with the r in previous equations representing the tank radius. Figure 10 shows how 4 varies with ~ and the frequency ratio. Also shown in Figure 10 is plot of K, the magnification factor. This factor enters into the solution equation and is evaluated as follows 1 K = (l-r ) + (2~r) where r is the frequency ratio. To apply the solution equation the constants K and 4 are first evaluated, then the constants of integration Al and A2 are determined based on the initial conditions. Then these values are used in the solution equation itself. 4. COMPUTER SIMULATION Now follows a discussion of the computer program developed to simulate the machine's operation. Figures 11(a) and (b) show a simplified flow diagram for the main program, a listing of which begins on page 36. Nte the diagrams shown in Figure 11(b) are connected through points and O). The program begins by making the necessary declarations and reading e user's input. Certain initializations are made and ISTART is set to 1. Additional input can be read if either INOUT or ITANK of the first input are set equal to anything but 1. Briefly, INOUT allows the user to control the program's output. It is possible to blank out certain portions of the output with this control. ITANK can be used if it is desired to start the run with a certain amount of mass in 23

0. 0. 0 1 2 3 Frequency ratio r 180o I 3 o. —-- — I -- - 0 12 30.7 45 0 I 2 3 Frequency ratio r Figure 10. Effect of damping on phase angle and magnification factor. 24

STARTSVl,INPfr INITIALIZATI ON =1W[} F I START = 1 ITANK=l (iffi) —-*|TA<TLEFT-TF T-" 0l) T F <im~sT LEFT(^o) T INITIALIZATION I._ _ SET=2 MEINEMTADDITIONAL ^ff ^*^, * 2 REVOLUTIONS INPU T 0'= CA LL TANK ISET T ~ FDW 5DA <^SET-1> --- ~ D Mt-DWAVE ~ ____ S - I Err IUITIAI TIME SPAN = SPEEDSEINTADERf CHANGE CODTOS+ ROTATION 10 I CO>NSTANTS j TAsTLEFT-TF SET CONSTANTS (l)- r'WD, FACT, PH I CON CA LL CA LC V1 I I ISTA RT _ _ TF=TO+TA CA LL 1 T TART I OLJTP1 DX<O CHANGE: CA LL IxIaX 00l.0 1^ TA NK MASS OUTP1 2 TOTAL MAS PROCEDURE FOR X>A."S. M^ TX NEGATIVE IS ^DX^a 10 SIMILAR. NO "^X'aOIOO(K ^y MASS CHANGES. TIME INCREMENT F = 300 ROTATION SET INITIAL CA LL CHANGE: T."^ CONDITIONS O — TANK MASS -(l + L=l TOTAL MASS FIND TA C104 } — FOR -(l0l) OX= (a) Figure 11. Main program flow diagram.

FIND TA ^ -^ C^O^^+ FFOR -- CA LL 0 00.02~^. X=AMP CA LC.^^ ~CHA NGE: T <T.FTE F TANK MASS I TOTAL MASS TF=TO+TA SET | C. SET 1. C. CA LL n A-<109)T |TIME SPAN TATLEFT-TF 110 _ —= — 20 ROTATION N CALL CAL T CX GF CHANG E: T ( ^112r< *,) ^X> XMA ) TANK MASS X>AMP l TOTAL MASS L CAFLL (115n \ OUTP2 MAX-2 IXMAX=X DX TA<TLEFT-TF TFMAXITF FIND TA TO=TF TF"MAX=TF FOR'~ —-T —-'~~~~~~~ m_____X:AMP CHA NGE: (IloG CALL._ TANK MASS I CALL OUTP1 ITOTAL MASSF CALC /TF=TO+TA (b) Figure 11. (Concluded)

the ballast tanks. Normally, however, INOUT and ITANK will be equal to 1 and this additional input will be skipped. The first time increment is taken as two rotations of the tank. This is to get past that portion of the time response where the transient is significant and into the smooth steady state response (see Figure 9(b)). Next the constants required by the solution equation are calculated. FACT and PHI are K and P described above. CALC is the subroutine called which takes the initial conditions and the constants just evaluated and determines the system response at the end of the given time increments and returns this information to the main program. The final time, TF, is set equal to the initial time, TO, plus the time increment, TA. If the program has just started and ISTART = 1, it outputs the information through subroutine OUTP1. ISTART is now set to 2 indicating the initial large time increment has been taken. At ) it is determined whether or not to increment the angular speed of the tank. If the speed, W, is already greater than the top speed desired, TOPW, the speed change routine is skipped. Also ISET must be equal to 1 for the speed to increase. ISET will be 1 only if the displacement of the tank on the previous half cycle was less than the desired amplitude, AMP. TANK is a subroutine which determines which ballast tank is adjacent to the flow nozzle and how much longer it will be there. The program sets the initial conditions and constants, and takes a time increment equal to the number of degrees of rotation as specified by input variable DEGREE. It makes sure, however, that this increment will not overlap two of the tanks. It does this by making sure TA < TLEFT - TF where TLEFT is the time at which the end of the current ballast tank will be reached, and is determined in the subroutine TANK. It then returns to 00) and the same course traveled before. This time, however, ISTART f 1 so it goes through the three conditionals. Normally it will pass through all three, update the mass in the ballast tanks, output the information, reset the initial conditions, and arrive back at. This process continues until the tank displacement becoms negative at which point the third conditional is n satisfied, and at e0 it finds the time at which X = 0 and returns to (10] or a final iteration. When X = 0 it switches out of the main loop and at ) updats the masses in the ballast tanks, outputs the information, and proceeds to 200) where the procedure for negative displacement begins. As indicated on th'diagram, this procedure is similar except the ballast tank masses do not change. The other condition, where it switches out of >e main loop, is if X > AMP. AMP is the desired displacement amplitude. o) connects to Figure 11(b). This portion of the program first accurately establishes the time at which X = AMP. It updates the masses in the ballast tanks and proceeds taking a 2~ of rotation time increment. The following steps are similar to those described above. This continues until X < AMP at which point it locates the X = AMP time, calls on CALC, updates the masses, and outputs the information prior to returning to -. This portion of the program is also set up to find the point at which X is the maximum (XMAX) and output this information. 27

A. SUBROUTINE CALC The flow diagram for this subroutine is shown in Figure 12(a). The calling variables for CALC are TA, X, DX, AREA. DX is the velocity X and AREA is the area under the displacement curve over the time increment TA. The necessary initial conditions and constants are obtained through the COMMON/VAR statement. The angular position of the tank is determined at To. The transient trig terms are evaluated along with the constants of integration Al and A2. Then after the steady state trig terms and EXP the exponential decay term are evaluated the solution equation is applied to find X, DX. AREA is determined by trapezoidal integration. This information is then returned to the main program. B. SUBROUTINE TANK Figure 12(b) shows a flow diagram for TANK whose calling variables are TF, TSAVE, W, ADJ, MTANK, TLEFT. These variables are defined as follows: TF - current time in the simulation run TSAVE - time at which present speed was reached W - present angular speed MTANK - a return variable indicating the current ballast tank TLEFT - time at which flow stream will start filling the next ballast tank First the position of the tank in fractions of a revolution- from the normal position (see Figure 4) is determined as REM. Then according to whether REM < 1/3 or 1/3 < REM < 2/3 or REM > 2/3 the ballast tank is determined as #1, #2, or #3. MTANK is set to this integer value, TLEFT is determined and the information is returned to the main program. C. SUBROUTINE FIND Figure 13 shows the flow diagram for FIND. This subroutine is called to interpolate a given table of values for X and TF to locate the time at which X has a certain value. The calling variables are: 28

SUB. SUB. CALC TANK (TA, X, DX, (TF, TSAVE, W, ADJ, AREA) MTANK, TLEFT) OBTAIN I.C. DETERMINE + CONSTANT THRU ANGULAR POSITION COMMON/VAR OF TANK, REM. I —— nCTCQtA.MC —— I^ rMTANK =1 DETERM I NE =T ANGULAR POS F T LEFT =T I M OF TA NK, BETAREMAO. 333 REMA INING FOR A T To FLOW INTO I AT l TQ' TANK 1 EVALUATE 1 TRANS IENT TRIG TERMS MTANK=2 /__ i ___> s F TLEFT=T IME ^EVA7LUATE 1REM>O. 666 REMA IN ING FOR TRANS lIENT FLOW CONSTANTS TANK 2 Al, A2 EVALUATE STEADY STATE TRIG TERMS MTANK=3 + EXP TLEFT=TIME REMA INING FOR FLOW INTO TANK 3 EVA LUATE X DX AREA 10 RETURN (a) (b) Figure 12. Subroutine flow diagrams. 29

SUB. FIND (XFIND, Z, F. L, Y) R1=IZ(L)-XFINDI RZ-=IZ(L-Z(L- I DT= I(F( L)-F(L-1) I1 R1, Y=F(L)+4- DT R2 RETURN Figure 13. Subroutine flow diagram.

XFIND - the value of displacement desired Z - table of displacement values F - corresponding table of times L - subscript of the last element of the above tables. The desired point is between elements L and L-l. Y - the returned value of time corresponding to XFIND The subroutine does a linear interpolation to find Y according to Y = F(L) + Z(L)- XFINDL) IF(L) - F(L-1)1 JZ(L) - Z(L-l) Y is then returned to the main program. D. SUBROUTINE OUTP1 This subroutine has a flow diagram as shown in Figure 14(a). It is called with the variables X, DX, TF, and TOTM where the latter is the total system weight. Output variables are obtained through the COMMON/VAR statement. Output control variables are obtained through the COMMON/OUT statement. As mentioned before these variables allow the user to blank out certain portions of the OUTPUT. If = 0 as preset in the main program, this output control is by-passed to (300). If Tl f 0, the next two conditionals let pass only times which are no tetween Tl and T2 or T5 and T4. Output between these times is not printed out as the program skips to 3). When it is desired to print the output, however, the output variables are TF, X, DX, BT(1), BT(2), BT(3), the masses in the three ballast tanks, TOTM, the total system mass, RPM, machine speed, and PHI, the phase angle between force and displacement. E. SUBROUTINE OUTP2 This subroutine shown by the flow diagram in Figure 14(b) is identical to OUTP1 except only TF, X, DX, RPM, and PHI, are output variables. This output subroutine is called during negative displacement when no mass changes occur, and during positive displacement when displacement is at a maximum and is above the AMP level. This helps to identify this point in the output listing. 31

SU S UB.U OUTP1 OUT P2 (X, DX, TF, (X DX TF TOTM) OBTAIN OUTPUT OBTA IN OUTPUT ARIABLES THRU VARIABLES THRU COMMON/VAR COMMON/VAR BTAIN OUTOBTAIN OUTPUTP CONTROL THRU CONTROL THRU COMMON/OUT COMMON/OUT T T T-=0 T1=0 I OP U UT 1R PM, P I /TFX, DX, BT(1),2 TF, X, DX, /BT(2), BT(3), TOTM,/ RPM, PHI RPM, PHI I ETURN) ETURN (a) (b) Figure 14. Subroutine flow diagram. 32

The following is a list of input and output variables with a description of their meaning. F. INPUT VARIABLES RPM - starting speed in rpm OFFBAL - weight of off-balance load in lb G4 - position of OFFBAL load measured counterclockwise in radians from ballast tank 2 TFIN - final time of run in sec DEGREE - amount of rotation in one time increment during the positive X half cycle. Should be less than 15~ ZETA - damping factor of the system INOUT - output control variable. If INOUT = 1, entire time response is output. If INOUT = any other integer, portions of output can be skipped using NAMELIST/OUTPUT FLOW - maximum flow rate in lb/sec TOTWT - total system weight in lb SPRING - spring constant in lb/in. RPMPS - acceleration in rpm/sec TOPRPM - top speed in rpm AMP - max amplitude of X in in. when X| is greater than AMP, m = FLOW and/or acceleration = 0 ITANK - input control variable. If ITANK is set equal to any integer but 1, starting mass in the ballast tanks can be input through "TANKIN" G. OUTPUT CONTROL If INOUT is set not equal to 1, Tl and T2, and T3 and T4 may be read in "OUTPUT." Output between Tl and T2 sec, and T5 and T4 sec will not be printed. 55

H. OUTPUT COLUMNS TIME - in sec (TF) X - displacement of the tank center in in. DX - velocity of the tank center in in./sec TANK 1 mass in tan bsecin. (BT( TANK 1 - mass in tank 1, lb-sec /in. (BT(l)) 2 TANK 3 - mass in tank 2, lb-sec /in. (BT(2)) 2 TANK 5 - mass in tank 5, lb-sec /in. (BT(5)) 2 TOT MASS - total mass of the system, lb-sec /in. (TOTM) RPM - speed in rpm PHI - phase angle between force and displacement, in radians The following is a list of other important variables in the program listing beginning on page 36, taken in order of first occurrence. RIN - input variable DEGREE in radians TSAVE - reset to the current time whenever the speed is changed ADJ - the angular position of the tank reset whenever the speed is changed UNBAL - OFFBAL converted to mass (lb-sec /in.') R - tank radius (in.) 2 DWSAVE - acceleration (rad/sec ) TOPW - TOPRPM converted to (rad/sec) ISET - memory interger controlling acceleration. Normally ISET = 1 and acceleration can take place. ISET = 2 if AMP was reached on prior half cycle and acceleration is not permitted. (Physically air clutch has been bled) ISTART - if set equal to 1 indicates run is beginning 54

DW - angular acceleration (rad/sec2) W - angular speed (rad/sec) TIN - time increment (sec) WN - natural frequency of system (rad/sec) RATIO - frequency ratio MTANK - number of current ballast tank TA - time increment to be taken (sec) TO - initial time at beginning of increment (sec) DXO - initial velocity (in./sec) XO - initial displacement (in.) WIN - speed increment (rad/sec) WD - damped natural frequency (rad/sec) FACT - magnification factor, K PHI - phase angle between force and displacement, 4 AREA - integral of the displacement curve over the time increment BT(MTANK) - mass in the ballast tank indicated by MTANK XFIND - argument of the value to be found by subroutine FIND RIN2 - angular increment taken during negative displacement (rad) I. OUTPUT An example of the output obtained is on page 44. The two types of format help to distinguish the positive and negative portions of the time response cycle. Note that there are many more interations taken during the positive displacement portion than there are during the negative. 35

PROGRAM IISTING $S I 5CX3 T=6C P=5C **LAST SIGNON WAS: 98:19.12 r4-16-71 USER "5nX3" SIGNED ON AT 12:27.14 CN 64-17-71 $LIST SPIN 1 NAMELIST/ ATA/RPMOFFBAL,G4,TFIN,DEGREE,ZETA, INOUT,FLOW, 2 1 TOT T,SPRING RPMPS,TCPRPPAMTP,IT ANK 3 REAC(5,DATA) 4 CIE\S ION Z(50),F(50 ),ET(3),BBT(3) 5 COMMCN/VAR/TO,TSAVEWPHI G4, XO CON,BTUNBAL,ZETA,WN,DXO 6 1 WD,ADJ,DW 7 COMMCN /OUT/T 1,T2,T3,T4 8 9 EQUIVALENCE(BBT ( 1 ),BT(1) ),(BBT( 2)BT(2)), ( BBT(3),BT( 3)) 10 RI N=DEGREE*6.28319/360(.0 11 WRITE(6,2) 12 2 FORMAT(' TIME',6X,1 X',SX,' DX',3X,' TANK l',3X,' TANK 2',3X, 13 1' TANK 3',3X,' TCT MASS RPM',5X,' PHI') 14 TSAVE=-. 15 CATA BBT,X,DXTF/(e.,0. Ti., 0,!,/ 16 T1=r. 17 FL=FL^/( AMP*3 E6.) 18 ADJ=n.9 19 UNEA L=J-FF AL/3 86. 20 TOTM=TOTWT/386. 21 P=12. 22 DWSA VE=R PMPS*6. 28319/6C. f 2 3 TOPW =T OPRPM*6. 2831 9/ 0.0 24 L=1 25 ISET=2 20 ITEST1l= 27 ITES12=1 28 ISTART=l 29 Dw O.0 30 IF(INOUT.EQ.1) GO TC 5 31 NAMEL IST/UUTPUT/ T1,T2,T3,T4 32 READ(5,OUTPUT) 33 5 CONT INUE 33.1 IF( ITANK.EQ.1) GU TO 6 33.2 N AELIST/TANKI /BT 33.3 READ(5,TANKIN) 33.4 6 CONTINUE 34 h = (RPi/6 3. )*6.2 8319 35 TIN=RIN/W 36 WN=S RT(SPRING/TOTM) 37 RATIC=W/WN 38 MTANK=1 39 TA= (2. *6.t) }/RPM 40 TO=TF 41 DXO=DX 42 XO-X 43 GO TC 1005 44 10 IF(TF.GE.TFIN) GO TO 3(00 45 IF(W.GE.TOPW) GO TO 11 46 IF(ISET.EQ. } GO TO 12 47 11 DW=O.e 48 GO TC 100 49 12 CW=Di SAVE 50 i IN=CW*TA 51 WNEW=W+WIN 52 REV=('TF-TSAVE ) *W/6.28319 36

53 IREV=REV+ACJ 54 REM= REV+ADJ-IREV 55 ADJNEW=REM 56 TIN RIN/WNEW 57 L=l 58 TSAVE=TF 59 W=WNEW 60 ADJ=ADJNEW 61 100 N=SQRT ( SPR I NG/TCT) 62 CALL TANK(TF,TSAVE,W, t'DJ,MTANK TLEF T) 63 RATIC=W/WN 64 TU=TF 65 CXC=CX 66 XO=X 67 TA=TIN 68 T ILEFT=TLEFT-TF 69 IF(TA.GT.TILEFT) TA=TILEFT 70 1005 D=(SQRT( 1.0-ZETA**2 ))* N 71 FACT=1./(SQRT( (1 O.-RAT IO *2) **2+( 2*Z ETA*RATIO } **2) ) 72 PHI=ATAN((2.0*ZETA*RATIC)/(1.O-RAT I**2 ) 73 CON= (RW*W*FACT )/SPR ING 74 101 CONTINUE 75 CALL CALC(TAX,DX,AREA) 76 TF=TC+TA 77 IF( ISTART.EQ. ) GO TO 1016 78 IF(DX.LT.~.i.AND.ABS(X).LE,.~.0001) GO TO 105 79 IF(X.GT.AMP.AND.ODXGT.. *0005) GO TO 106 80 IF(X GT.~. OC.AND.L.E. 1) GO TO 102 81 GO TO 104 82 1016 CALL OUTP1(X,DX,TFTCTM) 83 ISTART=2 84 GO TO 10 85 1012 BT ( MTANK)=BT (MTANK)+FL 4AREA 86 TOTN=TOTM+FL*AREA 87 CALL OUTP1(X,DX,TF,TGTM) 88 L=1 89 XO=X 90 CXO=CX 91 IF(ISET.EQ.1) GO TO ln 92 IF( TF.LT.TLEFT,AND. ISET.EQ.2) GO TO 103 93 GU TO 10 94 103 TILEFT=TLEFT-TF 95 IF( TA.GT.TILEFT) TA=TILEFT 96 TC = TF 97 GO TC 101 98 104 Z(L)=X 99 F(L) = TF 100 TA= 1A-. 01745/W 101 IF(X-.GE.0i.O) GO TO 1C45 102 L=L+1 103 CALL CALC(TA,X,DXAREA) 104 TF=TO+TA 105 GO TC 104 106 1045 XFIND=q0. 107 CALL FIND(XFIND, Z FL,Y) 108 TA=Y-TO 1C9 L=1 110 GO TC 101 111 105 BT(MIANK) =BT( MTANK)+FL*AREA 112 TOTM=TOTM+FL* AREA 37

113 CALL OUTP1(XDXTF,F,TCTt) 114 IF(IIEST1.EQ.1) ISET=1 115 ITEST2=1 116 GO TC 200 117 106 MAX=1 118 XMAX =.O 119 IF(X. LE(AMP+0.0002)) GO TO 1085 120 107 Z(L)=X 121 F(L)=TF 122 TA=TA-0 01745/W 123 IF(X.LE.AMP) GO TO 1C8 124 CALL CALC(TA,XOX,AREA) 125 TF=TO+TA 126 L=L+1 127 GO TC 107 128 108 XFIND=AMP 129 CALL FIND(XFIND, ZF,L,Y) 130 TA=Y-TO 131 = 1 132 CALL CALC(TAX,DX,AREA) 133 1085 BT(MTANK) BT (MTANK)+FL*AREA 134 TOTM=TOTM+FL*AREA 135 TF=TC+TA 136 109 TA=0.O349/W 137 XO=X 138 CXO=CX 139 TO=TF 140 L l 141 CALL TANK(TF,TSAVEW,ADJ,MTANK,TLEFT) 142 TILE FT=TLEFT-TF 143 IF( TA.GT.T ILEFT TA=TILEFT 144 110 CALL CALC(TA,X,DX,AREA) 145 TF=TC+TA 146 AREA=AMP*TA 147 IF(X.GE.AMP) GO TO 111 148 GO TO i15 149 111 BT(M ANK)=BT( MTANK)+FL*AREA 150 TOTM=TOTM+FL*AREA 151 IF(X.GT.XMAX) GO TO 114 152 IF(X.LE.XMAX.AND. MAX.GT.1) GO TO 112 153 CALL OUTP2(XMAX,DXMAX,TFMAX) 154 MAX=2 155 112 L=1 156 XO=X 157 DXO=CX 158 IF(TFLT.TLEFT) GO TC 113 159 GO TC 109 160 113 TI LEFT=TLEFT-TF 161 IF(TA.GT.TILEFT) TA —TILEFT 162 TO =TF 163 GO TC 110 164 114 XMAX=X 165 CXMAX=DX 166 TFMAX=TF 167 GO TO 112 168 115 CONTINUE 169 116 Z(L)=X 170 F(L)=TF 171 TA=T A-O.01745/W 172 IF(X.GE.AMP) GO TO 117 3&

173 CALL CALC(TA,X,CX,AREA) 174 TF=TC+T A 175 L=L+ 1 176 GO TC -116 177 117 XFINC=AMP 178 CALL FIND(XFIND,ZF,L,Y) 179 TA=Y-TO 180 ISET=2 181 ITEST1=2 182 L=l 183 CALL CALC(TAX,CX,AREA) 184 TF=TO+TA 185 AREA=AMP*TA 186 BT(MTANK) =T (MTANK)+FL*AREA 187 TOTM=TOTM+FL*AREA 188 CALL OUTP (X,DX F,TF,TCT) 189 GO TO 10 190 200 RIN2=0.5236 191 ILEFT= TF+(3.14159/W) 192 210 IF(TF.GT.TFIN) GO TO 3C0 193 TA=RIN2/W 194 IF(W.GE.TLPW) GC TO 211 195 IF(ISET.EQ.1) GO TO 212 196 211 CW=.0 197 GO TO 220 198 212 OW=D SAVE 199 I N=CW*TA 200 REV=(TF-TSAVE) */6.28319 231 IREV=REV+ADJ 202 REM= REV +ADJ-IREV 203 ADJ=REM 234 W=W+WIN 205 TA=RIN2/W 206 L=l 207 TSAVE=TF 208 220 WN=SCRT (SPRING/TGTM) 209 RATIC=W/WN 210 CXO=CX 211 XG=X 212 TO=7F 213 221 WD=(SQRT(l. -Z ETA**2)) 4N 214 FACT=1.0/(SQRT((.n-RATIC**2)**2+(2*ZETA*RATIO)**2)) 215 PHI=ATAN((2.0*ZETA*RATIC)/(1.0-RATIO**2 ) 216 CON= (R*W*W*FACT )/SPRING 217 222 CALL CALC(TA,X,CX,AREi) 218 TF=TO+TA 219 IF(DX.GT.n.,O.AND. ABS(X).LE..000 N1) GO TO 236 220 IF(X.LT.-AMP.AND. CX.GT.O.0005) GO TO 227 221 IF(X.LT.O.O.AND.L.EQ.1) GO TO 223 222 GO TC 225 223 223 CALL OUTP2(X,DX,TF) 224 TO=TF 225 XO=X 226 DXC=DX 227 IF((TF+.5236/W).GT.TLEFT) GO TO 224 228 IF(ISET.EQ.1) GO TO 210 229 GO TC 222 230 224 TA=3.08727/W 231 GO TC 222 232 225 Z(L)=X 39

233 F(L)=TF 234 TA=TA-."^1745/W 235 IF(X.LL.(.O ) GO TO 226 236 L=L+i 237 CALL CALC(TA,X,CX,AREA) 238 TF=TO+TA 239 GO TO 225 240 226 XFIN C=". 241 CALL FIND(XFIND,Z,FL,Y) 242 TA=Y-TU 243 L=l 244 GO TO 222 245 227 MIN=1 246 XMIN=r=.0 247 IF( X.GE.-(AMP+0.f.C)02) 1 GC TO 230 248 228 Z(L)=X 249 F(L)=TF 250 TA=TA —.01745/i 25! IF(X.GF.-AMP) GO TO 229 252 CALL CALC (TA,X,DX, PREA 253 TF=TO+TA 254 L=L+l 255 GO TC 228 256 229 XFIND=-AMP 257 CALL FINO(XF INOZ,F,L,Y) 258 TA=Y-TU 259 L=l 26 ) CALL CALC(TAX, CX, AREA ) 261 TF=TC+TA 262 230 1A =:.' 349/ 263 XO=X 264 DXO=DX 265 TO=TF 266 231 CALL CALC(TA,X,DX,AREA) 267 TF=TC+TA 268 IF( X.GT.-AMP) GC TC 234 269 232 ClONT INUE 27 IF(X.LT.XMIN) GC TC 2?3 271 IF(X.GF.XIMIN.ANO. MIN.GT.1 GO TO 230 272 CALL UlUTP2 (XM I N, DX i I N, TFM IN) 273 MI=2 274 GOJ TC 23C 275 233 X IN=X 276 D)X =OX 277 TFMI N=TF 278 GO TC 23: 279 234 Z(L)=X 2W'n F(L)=TF 281 TA= T -..!' 1745/ 282 IF( X.LLE.-AIP) GC TC 235 283 CALL C!ALC(TA,X, UX, A EA) 2 4 L=L+1 2b5 GO TO 2J4 286 235 XF I i4 C= - AMP 287 CALL FINf)(XFING, Z, FL,Y) 288 TA=Y-Tr 289 ISET=2 2 9I T E S T 2 = 291 L=l 292 CALL CALC(TA,X,CX, AiFA ) 40

293 CALL OUTP2(X,DX,TF) 294 GU TO 210 295 236 CALL OUTP2(X,DX,TF) 296 IF( IIEST2.EQ.1) ISET=1 297 ITEST.=1 298 L=1 299 GO TC 10 300 300 CONTINUE 301 END END OF FILE 41

$LIST SUB 1 SU RCUT IE CALC(TAXtDXAREA) 2 LIMENSION BT(3) 3 CCMMCN/VAR/TO,,SAVE,W,PFIG4, XOCON,BT,UNBALZETAWN,DX0, 4 1 hD,ADJ,DW 5 RAC= (TO-TSAVE )*W 6 REV=RAD/6.28319 7 IREV=RE V+ADJ 8 REM=REV+ADJ-IRE" 9 BET A=REM*o 28319 10' THETA1=BETA-2.618-PH I 11 THETA2=BETA+1.57C8-PHI 12 THET A3=BtETA-O.52 36-P-F I 13 THE TA4=BE TA+1. 57T08-G4 —PH 14 S1=SIN(THETA1 ) 15 S2=SIN(THET A2 ) 16 S3=SIN(THETA3) 1 7 S4= S IN(THJT A4) 18 CLN2 =( CON*DW ) / ( b**2 ) 19 C1=COS( THETA1) 20 C2=CCS THFTA2) 21 C3=CCS(THETA3) 22 C4=CS( THETA4 ) 23 A1=XC-CON* ( BT ( 1 )*S 1 +T (2 )*S2+ BT ( 3) *S3+UNBAL*S4) +CN2* 24 1 (BT( 1)*C1+BT( 2)*C2+T( 3)*C3+UNBAL*C4) 25 A2=( A1*ZETA*WN+DOXU-CCN*W*(BT( )*C1+BT( 2)*C2+BT ( 3)*C3+UNBAL*C4 ) 26 1 /W0- CUN2*W* (T( 1)*S1+tT (2)*S2+BT(3 )*S3+UNBAL*S4) ) /WD 27 WT=W*TA 28 SW1=SIN(WT+THETAl) 29 S^2=SIN(WT+THETA2 ) 30( SW'3= SI N( T+ T HE TA 3) 31 S4=SI N (WT+THETA4) 32 CW1=COS (T+THETA1) 33 CW2=COS(WT+THETA2) 34 CW3=COS(WT+THETA3) 35 CW4=COS(( T+THETA4) 36 SWE=S I4( WD*TA ) 3 7 C D=COS (WD*T A) 38 EXP=2.71828**( ZETA*N*TA) 3 9 STER =BT ( 1 ) SW 1+ BT ( 2 ) SW 2+BT( 3 )*SW3+UNBAL*SW4 40 CTER =PT(1 )*CW1+BT(2)*CW2+BT (3)*CW3+UNBAL*C W4 41 X=( A 1*C D+A 2*SWD )/E XP+CON*STE RM-CN 2* CTERM 42 CX=(-ZETA*WN*(Al1*CWO+A2*SWD )+WU*(-Al*SWD+A2CWO ) /EXP+CON*W* 43 1 CTERM+CON2*W*STERM 44 AREA=.5*(XO+X ) *TA 45 RETUFN 46 END 47 SUBROUTINE TANK(TF,TSAVEt,,ADJ,MrANK,TLEFT) 48 49 REV= ((TF-TSAVE.)*;W)/6.28319 50 IREV=REV+ADJ 51 REM=REV+AJ- IPEV 52 IF(REM.GE.0.999) REM=..() 53 IF(REM.GE.O.333) GO TO 11 54 PTANK = 1 55 TLEFT=TF+( <.. 3333-R EM)*6.28319 / 56 GO TC 1r' 57 11 IF(RE:.GE.O'.666) GO TO 12 58 MTANK = 2 42

59 TLLFT =T F+(I.6667 -REM ) *6. 28319/W 60 GO TC 1)i 61 12 MTANK = 3 62 TLFFT=TF+(1.0 -REM)*6.28319/W 63 100 CONTINUE 64 RETU RN 55 END 66 SUBROUTINE OUTP1(XtDX, IFTCTM) 67 - DIMENS ION B (3) 68 CCMMCN/V AR/TC, TSAVE,, Pf I,G4 XO, CON, BT UNBAL, ZET AWNt DX0, 69 1 WD,ADJ,DW 70 CCMMCN/OUT/T1, T2,T3 T4 71 IF(Tl.EQ.O.O) GO TO 3(0 72 IF(TF.GT.T1.AND. TF.LT.T2) GO TO 450 73 IF(TF.GT.T3.ANDO TF.LT.T4) GO TO 450 74 3C 0 CONTINUE 75 RPM=W*60.. 1/6.28319 76 hRITE(6,4'r) TFX,OX,8TI),T (2), BT(3),TOTMRPM,PHI 77 4r" FORMAT(F1'0.4,2F1Cn.5,4F1C.7,F8.3,F8.4) 78 4'5 CONT \NUE 79 RETURN 80) END 81 SUBROUTINE OUTP2(XtDXTF) 82 CIME\SION ET(3) 83 COMMON/VAR/TO,TSAVEWPHI,G4,XUCONBTUNBALZETA,WNDXO, 84 1 WD,ADJDW 85 COMMCN/OUT/T 1 T2,T3 T4 86 IF Tl.EQ.0.CI GO TO 3#-.0 87 IF(TF.GT.Tl.AND. TF.LT.T2) GO TO 450 88 IF(TF.GT.T3.AND. TF.LT.T4) GO TO 450 89 3(-0 CONTINUE 90 RPM=W*60.C/6 28319 91 WRITE(6,410) TF,X,DX,RPFPHI 92 410 FORM AT(Flr0.4,2F10.5,4i0X,F8.3,F8.4) 93 450 CONTINUE 94 RETURN 95 END 96 SUBRCUT INE FIN ( XFINC ZFLY ) 97 DIMEtNSION Z( 5O) tF(5) 98 R1=AES(Z(L )-XF IND) 99 LL=L-1 100 R2=AES(Z(L)-Z(LL ) 101 CT=ABS ( F(L )-F(LL L ) 102 Y=F( )+(Rl/R2) *DT 103 RETURN 1)4 END ENO OF FILE 45

OUTPUT LISTING $RUN SPINOBJ+SUBOBJ 5=*S3URCE* 6=*SINK* EXECUTION BEGINS TIME X DX TANK 1 TANK 2 TANK 3 TOT MASS RPM 0.8030 0.00033 0.15379 0.0 0.0 0. 0 0,5181347 150.000 0.8133 0.00236 0.14945 0.0000031 0.0 0., 0.5181378 150.000 0.8267 0,0042~ 0.13833 0.0000106 0.0 O,0 0,5181453 150,000 0,8430 0.00602 0.12086 0.0000224 0.0 0.0 0,5181571 150,000 0.8533 0.00748 0.09807 0.0000378 0.0 0.0 0.5181724 150.000 0.8667 0.00862 0.07121 0.0000561 0.0 0.0 0.5181907 150,000 0.8800 0,00937 0.04153 0.0000766 0.0 0.0 0.5182112 150.000 0.8933 0.00971 0.01019 0.0000983 0.0 0,3 0.5182329 150,000 0.9067 0.00963 -0.02155 0.0001204 0.0 0.0 0.5182549 150.000 0.9200 0.00914 -3.05229 0.0001418 0.0 0.0 0.5182762 150.000 0.9333 0.00825 -0.08057 0.000616 0.0 0.0 0.5192960 150.000 0.9457 0.00701 -0.10524 0.0001616 0.0000174 0.0 0.5183133 150.000 0.9600 0.00546 -0.12584 0.0001616 0.0000316 0.0 0.5183275 150,000 0.9733 0.00367 -0.14148 0.0001616 0.0000420 0.0 0.5183379`150.000 0.9867 0.00171 -0.15070 0.0001616 0.000481 0.0 0.5183440 150,000 0.9979 -0.00000 -0.15284 0.0001616 0.0000498 0.0 0.5133456 150.000 1.0311 -0.00475 -0.13396 150.637 1.0642 -0.00857 -0.08237 151.271 1.0971 -3.00995 0.00008 151.902 1.1298 -0.00878 0.07639 152.530 1.1625 -0.00516 0. 14192 153.155 1.1950 -0.00011 0.16360 153. 780 1.1957 -0.00000 0.16368 153.780 1.1971 0.00023 0.16379 0.0001616 0.0000498 0.0000000 0.5183456 153.792 1.2101 0.00235 0.16082 0.0001644 0.0000498 0.0000000 0.5183484 153.819 1.2230 0.00438 0.15057 0.0001719 0,0000498 0.0000000 0.5133558 154.067 1.2360 0,00622 0.13329 0.0001836 0.0000498 0.0000000 0.5183675 154.315 1.2489 0.00781 0.10984 0,0001991 0.0000498 0.0000000 0.5183830 154.562 1.2619 0.00905 0.08150 0.0002177 0.0000498 0.0000000 0.5184016 154.809 1.2747 0.00990 0.04961 0.0302386 0,0000498 0.0000000 0.5184225 155.056 1.2876 0.01032 0.01548 0.0002609 0.0000498 0.0000000 0.5134447 155.302 1.3005 0.01030 -0.01953 0.0002835 0.0000498 0.0000000 0.5184674 155.548 1.3133 0.00982 -0.05398 0.0303056 0.0000498 0.0000000 0.5184894 155.794 1.3261 0.00893 -0.08627 0.0003261 0.0000498 0.0003003 0.5185099 156.039 1.3389 3.00763 -0.11500 0.0003261 0.0000679 0.0000000 0.5185280 156.283 1.3517 0.00601 -0.13947 0.0003261 0.0000828 0.0000000 0.5185428 156.528 1.3645 0.00410 -0.15861 0.0003261 0.0000938 0.000000-0 0.5185538 156.772 1.3772 0.0019 -0.17084 0.0003261 0.0001004 0.o000000 0.5185604 157.015 1.3887 0.00000 -0.17494 0.0003261 0.0001024 0.0000000 0.5185623 157.258 1.4204 -0.00534 -0,15342 157.85o 1.4519 -0.00932 -0.09045 158.470 1.4833 -0.01085 -0.00079 159.073 1,5146 -0.00948 0.09060 159. 673 1.5458 -0.00552 0.15999 160. 271 1.5769 -0.000001 018734 160.867 1.5781 0.00022 3.18746 0.3003261 0.0001024 0.0000000 0.5185623 161.460 1.5905 0.00253 0,18435 0.0003290 0.0001024 0.0000000 0.5135652 161.484 1.6029 0.0.0475 0.17314 0.0003367 0.0001024 0.0000000 0.5185729 161.720 1.6152 0.00678 0.15352 0.0003489 0,0001024 0.0300000 0.5185850 161.956 1.6276 0.00851 0.12632 0.0003650 0.0001024 0.0330000 0.5186011 162.192 1.6399 0.00987 0.09336 0.0003843 0.0001024 0.0000000 0.5186204 162.427 1.6522 0.01079 0.05666 0.0004060 0.0001024 0.0 30000 0.5186421 162.662 1.6645 0.01125 0.01783 0.0004292 0.0001024 0.0000000 0.5186653 162.897 1.6767 0.01123 -0.02189 0.0004527 0.0001024.00300003 0.5186888 163.132 1,6890 0.01072 -0.06113 0.0004757 0.0001024 0,030000330 5137117 163.366 1.7012 0.00975 -09.09819 0.0004971 0.0001024 0.0000000 0.5187331 163.599 44

J. HOW TO USE THE PROGRAM The following are the commands required to do a typical run. It is assumed the program has been compiled and is in object file OBJ. For completeness the run will start with masses in the ballast tanks, and it will be desired to blank out some of the object. $RUN OBJ 5 = *SOURCE* 6 = *SINK* &DATA RPM = 150., OFFBAL = 5., G4 = 0., TFIN = 40., DEGREE = 10., ZETA = 0.1., INOUT = 3, FLOW = 0.695, TOTWT = 200., SPRING = 500., RPMPS = 2., TOPRPM = 450., AMP = 0.7, ITANK = 2 &END &OUTPUT Tl = 15., T2 = 30., &END &TANKIN BT = 0.002, 0., 0.004 &END $ENDFILE Note the input is name listed and starts in column 2 on the card. If INOUT = 1, this would mean no &OUTPUT card. If ITANK = 1, no &TANKIN card would be read. K. RESULTS AND CONCLUSIONS The results of two runs have been reduced to amplitude and speed vs. time. These are shown in Figures 15 and 16. The two runs are identical except for the flow rate. Figure 15 has a flow rate comparable to the present washer-dryer and in Figure 16 this has been increased by 50%. As might be expected, the higher flow rate allows the machine to attain maximum spin speed faster. In interpreting the results of various simulation runs, the differences between the actual machine and the computer model must be kept in mind. The model does not incorporate the air clutch speed control system. A direct simulation of this system was considered outside the scope of this project. Therefore, it was replaced by an ON-OFF constant acceleration system whereby the tank is either accelerating at the specified input rate or is not accelerating at all. The criterion being whether or not the displacement on the prior half cycle was greater than AMP. As a result it would not be expected that the simulation would yield the same time-speed characteristic as the actual machine, but comparison of two simulations as in Figures 15 and 16 might be expected to parallel changes in machine hardware. As described above it is possible to vary certain parameters and investigate the effects of hardware changes on the machine. It is also possible to investigate other possible spinning schemes. One such scheme which was tried was to balance the load at low speed and then spin it at a speed above critical. This involves two separate runs. First the off-balance load was balanced at 150 rpm for 40 sec with 5 gal/min flow rate. The suspension in 45

400.015 300 Speed RPM,010 -200 /Unbalance= 5 lbs. Flow = 5 lbs..005- Flow = 0.695 Iblsec.= 5 gal./min. K = 4040 Ib/in. \ IS Total Weight= 200 Ibs. E c I i I I 10 20 30 40 TIME SECONDS Figure 15. Time to attain balanced spin-data from program output.

-400.015300 Speed /RPM.010 I -200 Unbalance = 5 lbs. Flow = 1.0425 Ib/sec..005- 7.5 gal./min. K= 4040 lb./in. 100 Total Weight =200 lbs. Cn E- a 10 20 30 40 TIME -SECONDS Figure 16. Time to attain balanced spin-at a different flow rate. Data from program output.

this scheme was softened to 500 lb/in. which gave a critical speed of 300 rpm assuming the same machine weight. Then with the final masses in the ballast tanks given as input, a run was made starting at 150 rpm accelerating at 50 rpm/sec to 600 rpm with no flow rate. The value for t was.1 in this run and at critical an amplitude of.5 in. was reached. The limitation of the program in simulating a soft suspension such as this is, of course, the fact that the simulation has 1 degree of freedom. In a soft suspension system the tank will be allowed vertical as well as horizontal deflection. But it would be possible to make changes in the program to take into account the effects of coupling between the coordinates. Instead of calling on the subroutine CALC, a vibration program such as VIB which handles multi-freedom systems could be called. It would return the horizontal deflection which would then be fed back into the balance control as is presently done. Thus it would be possible to incorporate into the present ballast tank simulation a multi-freedom model. REFERENCES - Ogata, Katsnhiko, Modern Control Engineering, 1970, Prentice Hall, Inc., Englewood Cliffs, N. J., pp. 663-695. - Carnahan, B., and J. 0., Wilkes, Digital Computing and Fortran IV 1968, The University of Michigan, Ann Arbor, Michigan. - Chace, Milton A., Dynamics of Machinery Systems-A Vector/Computer Oriented Approach, 1970, Prentice-Hall, Inc., Englewood Cliffs, N. J., Chapters 5 and 6. - Tse, Francis S., Ivan E. Morse, and R. T. Hinkle, Mechanical Vibrations, Allyn and Bacon, Inc., Boston, 1963. 48

A COMPUTER SIMULATION OF A VENTLESS CLOTHES DRYER USING A REFRIGERATION CYCLE Kenneth J. Timmer 49

1. OBJECTIVES The purpose of this project was to develop a digital computer program to simulate the behavior of a closed-loop or ventless clothes dryer. This system passes the dryer exhaust air over the evaporator and condenser coils of a refrigeration cycle to change the properties of the air in the feedback portion of the closed-loop system (Figure 1). The principal reasons for developing this computer simulation were to determine the feasibility and capability of such a clothes dryer system and to aid in the actual design of such a system. 2. RESULTS The computer simulation program consists of a main program which performs the principle calculations, a psychrometric subroutine which calculates the thermodynamic properties of the air-vapor mixture, a plotting subroutine, and a somewhat modified version of the dryer simulation subroutine package develope' by Whirlpool. The general results of the simulation indicate that such a clothes dryer system is feasible. The system performs the primary objective of drying clothes and does so in a reasonable amount of time. The system has the advantages of not requiring any vent and operating on normal 110-volt household current. The simulation results show that although the drying time required is about 36O/ greater than with the conventional electric dryer which requires a 220-volt source and uses a 5600-watt heater, the amount of power required is only 1/3 of that required with the conventional system. Figure 2 compares the power and drying time required of a 5600-watt heater system with the simulated close-loop refrigeration system. The figure also includes a 1200-watt heater system which is approximately the wattage required by the refrigeration system, and a system using a refrigeration system but in an open-loop configuration. A schematic diagram of an open-loop system is given in the Appendix. 3. BACKGROUND The drying of clothes in an electric clothes dryer is accomplished by passing warm dry air over the clothes. The water in the clothes is evaporated and carried out of the dryer by the air stream. In a conventional dryer this 51

REGION REGION C DRYER A i /^- 8 mXEXPANSION VALVE REGION CONDENSER EVAPORATOR DRAIN COMPRESSOR Figure 1. Schematic diagram of ventless dryer system. 52

140 3.36 134.7 3 120 2.69 100 2 - 80 W W0 55.8 1.12 45.0 40 36.0 0.90 20 5600 Watt 1200 Watt Refrigeration Refrigeration System System System System (Closed Loop) (Open Loop) Figure 2. Comparison chart of power and time required to dry a normal clothes load (9 lb dry weight, 850 initial moisture content). 55

air is exhausted to the outside and a fresh quantity of dry air is constantly being drawn into the dryer. In order to operate a clothes dryer using a ventless system, some method is needed to remove water from the air stream to keep the air going over the clothes dry enough to pick up moisture from the clothes. To condense water the air-vapor mixture must be cooled below the dew point of the mixture. This can be accomplished with a variety of types of heat exchangers. However, once the air-vapor mixture has been cooled below its dew point and water has been condensed out to reduce its specific humidity (ibm H20/lbm air), it must then be reheated to reduce its relative humidity and enable it to continue to absorb moisture from the clothes. This requires another heat exchanger or electrical heating element. If these two heat exchangers are operated separately, the energy obtained when the air-vapor mixture is cooled is merely dissipated and a significant amount of additional energy is required to heat the air. This brings us to the investigation of using a refrigeration cycle in a ventless clothes dryer. Not only does a refrigeration cycle have the necessary heat exchangers-the evaporator for cooling and condensing water and the condenser for heatingbut it absorbs the energy obtained in cooling the air at the evaporator and turns it to energy used in heating the air at the condenser. Therefore, it is possible to operate such a system with only a fraction of the power required with separate heat exchangers. 4. BASIC EQUATIONS It is possible to analytically represent the changes in the air-vapor mixture as it passes over the evaporator and condenser coils of the refregeration cycle with energy balance equations once the amount of heat transfer at each coil is known. The terms which are necessary for these equations are: QE - heat transfer at the evaporator (Btu/hr) QC- heat transfer at the condenser (Btu/hr) m - mass flow of air (lbm/hr) h - enthalpy of air (Btu/lbm) X - specific humidity of air (Ibm H20/lbm air) 54~~

h - enthalpy of vapor (Btu/lbm) h - enthalpy of liquid (Btu/lbm) If we denote regions A, B, and C as in Figure 1, the basic equations of the simulation are: At the evaporator h +w hB + ( h ()h (1) A A A B B A A B ( ENTHA + -QEVAP = ENTHB + AE ENTHB = ENTHA - QEVAP - AE (2) At the condenser h ~"B h + Q m h (3) B B VB C B C C o ENTHB + QCOND = ENTHC ENTHC = ENTHB + QCOND (4) In these two equations the energy quantity of the air and vapor in any particular region can be lumped together to form an ENTH quantity which refers to the energy of the air-vapor mixture. These energy terms can be found directly using a psychrometric chart. The AE term is very small relative to the other terms in equation (2) and is assumed to be a constant fraction of ENTHB (see Appendix). The mass flow of the air can be calculated from the given volumetric air flow in ft3/min and the density of the air vapor mixture which is also found from the psychrometric chart. This leaves QE and QC as the only terms to be determined. During much of the initial work on the program these two terms were given constant values based on what heat transfer quantities could reasonably be expected with a compressor which operates on 110 volts. The values used during this work were 11,000 Btu/hr heat transfer at the evaporator (QE), and 15,000 Btu/hr heat transfer at the condenser (Qc) using Freon-22 as the refrigerant. 55

5. ANALYSIS OF THE REFRIGERATION CYCLE The type of refrigeration cycle used in the simulation was a vaporcompression cycle using a refrigerant such as Freon-22 as the working fluid. In such a system, heat is transferred to the refrigerant in the evaporator, work is done on the refrigerant in the compressor, and heat is transferred from the refrigerant in the condenser. The temperature-entropy diagram for this system is shown in Figure 3. The amount of heat transferred from the refrigerant to the air in the condenser is the sum of the heat transferred to the refrigerant from the air in the evaporator and the work done on the refrigerant in the compressor. In the closed-loop dryer the same air stream passes over the evaporator and condenser. Therefore, since the heat transferred to the air is always greater than the heat transferred from the air, the temperature of the air stream will be constantly increasing. Because of this constant temperature buildup it is erroneous to assume constant heat transfer quantities. Increasing air temperatures cause increasing refrigerant temperatures and these in turn cause changes in the heat transfer quantities and amount of compressor work required. Published compressor curves prepared by pump manufacturers show the effect of evaporator temperature and high side, or condenser, pressure on the capacity or heat transfer of the evaporator and on the power required in the compressor (see Appendix for example). The high side pressure in turn is directly related to the temperature of the refrigerant in the condenser. These compressor curves can be reduced to equation form to be used in the computer simulation. To determine the effect of rising air temperatures on the properties of the refrigerant and thus be able to calculate heat transfer quantities, a series of constraints were used in the computer program. With the use of this part of the program the properties of the air and refrigerant are constantly checked and as the air temperatures increase the refrigerant temperatures are incrementally increased and new heat transfer quantities are calculated and used in the basic system equations. The constraints used are: 1. The temperature of the refrigerant in the evaporator increases as the enthalpy of the air passing over the evaporator increases. This relationship is found from McQuay refrigeration coil data. A designed air enthalpy value is chosen by the user of the simulation. When the enthalpy of the air increases above this chosen value, the evaporator temperature increases at a rate determined from the McQuay data.

TEMPERATURE ENTROPY Figure 3. Temperature-entropy diagram of refrigeration system. EXPANSION VALVE CONDENSER EVAPORATOR 3 2 CO MPRESSOR Figure 4. Schematic diagram of refrigeration system. 57

2. The temperature of the refrigerant in the condenser increases as the temperature of the air over the condenser increases. A designed temperature for air going into the condenser is chosen by the user. When the air temperature increases above this value, the condenser temperature also increases at a rate determined from the McQuay data. 3. The amount of heat transferred in the condenser is limited by the temperature of the refrigerant in the condenser. A maximum temperature of 5~ less than the condenser temperature is defined. The maximum amount of heat transferred is then that amount which raises the air temperature to the maximum just defined. Since this has the effect of increasing the amount of energy remaining in the refrigeration cycle, an incremental temperature difference chosen by the user is added to the refrigerant temperatures of both the condenser and evaporator. 4. As indicated above, the condenser temperature and high side pressure are directly related. The high side pressure is limited by the capability of the compressor to pump the refrigerant above a certain maximum value without burning out. This maximum pressure is chosen by the user based on knowledge of the refrigerant and compressor being used. When the pressure in the condenser reaches this maximum value the capacity or heat transfer quantities that the evaporator and condenser have at that time are assumed to remain constant for the remainder of the drying time. Possible suggestions to accomplish this steady state period include bleeding enough air from the dryer exhaust and replacing it with room air to reduce the energy content of the air the same amount as it is increased as it passes over the evaporator and condenser (see Appendix for this calculation). This series of constraints was based on the best information I could find in the time available. It is quite possible, however, that as a prototype such a system is tested better relationships between air and refrigerant properties will be found. These relationships can then be put into the simulation without destroying the basic logic of the program. 6. DISCUSSION OF MAIN PROGRAM The simulation of this ventless clothes dryer is a combination of the basic equations discussed and the four constraints necessary in determining the capacities of the evaporator and condenser of the refrigeration cycle. 58

The input data which must be read into the program are: (1) TDB = Ambient dry bulb temperature [~F] (2) TWB = Ambient wet bulb temperature [~F] (3) CLOAD = Dry weight of the clothes load [lb] (4) CMOIST = Initial moisture content of the clothes as a percentage of the dry weight [%] (5) CTEMP = Initial clothes load temperature [~F] (6) CFM = Volumetric flow of the air-vapor mixture in the closed-loop [ft3/min] (7) HSP = Initial high side pressure of the refrigerant [psig] (8) EVTP = Initial temperature of the refrigerant in the evaporator [~F] (9) HSPMAX = Maximum high side pressure of the refrigerant [psig] (10) HMAX = Value of the enthalpy of the air-vapor mixture going into the evaporator used in the refrigeration cycle constraints [Btu/lbm] (11) TBMAX = Value of the temperature of the air-vapor mixture going into the condenser used in the refrigeration cycle constraints [~F] (12) CONADD = Temperature added to the temperature of the refrigerant in the condenser as part of the constraints [~F] (13) EVADD = Temperature added to the temperature of the refrigerant in the evaporator as part of the constraints [~F] (14) IPLOT = Control integer governing the number of plots printed out: 0 = no plots 1 = one plot 2 = two plots 59

The meaning of symbols of other properties used in the program are: (1) MMOIST = Moisture content of the clothes [ibm water] (2) WATTS = Wattage of the heating element [watts] (3) TDBA, TDBB, TDBC = Dry bulb temperature of the air-vapor mixture in region A, B, and C, respectively [~F] (4) TWBA, TWBB, TWBC = Wet bulb temperature of the air-vapor mixture in region A, B, and C, respectively [~F] (5) HUSA, HUSB, HUSC = Specific humidity of the air-vapor mixture in region A, B, and C, respectively [Ibm H20/lbm dry air] (6) DENSA, DENSB, DENSC = Density of the air-vapor mixture in region A, B, and C, respectively, [lbm dry air/ft3] (7) ENTHA, ENTHB, ENTHC = Enthalpy of the air-vapor mixture in region A, B, and C, respectively [Btu/lbm dry air] (8) MDAA, MDAB, MDAC = Mass flow of the air-vapor mixture in region A, B, and C, respectively [Ibm dry air/hr] (9) DENSEV, DENSCO = Average density of the air-vapor mixture in the evaporator and condenser, respectively [Ibm dry air/ft3] (10) MDAEVA, MDACON = Average mass flow of the air-vapor mixture in the evaporator and condenser, respectively [Ibm dry air/hr] (11) CONTP = Temperature of the refrigerant in the condenser [~F] (12) QOUT = Heat transferred in the evaporator [Btu/hr] 60

(13) POWW = Power required in the compressor [watts] (14) POWB = Energy required in the compressor [Btu/hr] (15) QADD = Heat transferred in the condenser [Btu/hr] (16) QEVAP = Energy transferred in the evaporator [Btu/lbm dry air] (17) QCOND = Energy transferred in the condenser [Btu/lbm dry air] (18) TMAX = Maximum temperature of the air out of the condenser [~F] (19) HDIFF = Enthalpy difference [Btu/lbm dry air] (20) TADDE = Temperature to be added to the refrigerant in the evaporator [~F] (21) TADDC = Temperature to be added to the refrigerant in the condenser [~F] (22) CONMIN = Minimum temperature of the refrigerant in the condenser [~F] (23) EVTPIN, CONIN = Initial temperature of the refrigerant in the evaporator and condenser, respectively [OF] (24) TIME = Actual time the dryer has been operating [hr] (25) DRATE = Rate at which water is being removed from the clothes [Ibm H20/hr] (26) HOR = Variable on the horizontal axis of the printed plots (27) FUNC1, FUNC2 = Variables on the vertical axis of the first and second plots, respectively PROGRAM DISCUSSION Numbers correspond to numbers on the flow diagram which follows. (1) To begin the program the input data are read in. These values can be put on cards or kept in the memory file. 61

(2) The initial moisture content of the clothes is calculated from the input data. (3) The DRYER subroutine requires a value for WATTS, but since this simulation does not include a heating element the variable is set equal to zero. The air going into the evaporator is set equal to the ambient dry and wat bulb temperatures. (4) The PSYCKT subroutine is called to calculate the other properties of the air in region A and the mass flow of air is calculated. (5) Enough information has been calculated to call the DRYER subroutine to initialize the time and give initial values to the properties of the clothes. The DRYER subroutine does not perform any calculations to simulate drying of the clothes at this time. (6) Initial values of zero are given to all the control integers and initial values of the temperature of the refrigerant in the evaporator and condenser are calculated. The relationship of condenser temperature to high side pressure would differ with a different compressor or refrigerant. (7) The amount of heat transferred from the air-vapor mixture to the refrigerant in the evaporator is calculated: QOUT = 50(300 - HSP) +.281(800)(EVTP - 35) + 9200 This equation is one of a series of equations derived by the Air Conditioner Division of Whirlpool in Evansville, Indiana from actual compressor curves. This particular equation is for a Tecumseh B1613 compressor using Freon-22. (8) Based on the evaporator temperature and high side pressure, the power required at the compressor is calculated, using a relationship derived from this particular compressor curve. POWW = [4.35 +.o6(HSP - 260)](EVTP - 43.5) + 1.75(HSP - 260) + 1030 This equation is accurate when: EVTP > 43.5 ~F, and HSP > 260 psig 62

(9) The power at the compressor is limited to a maximum of 1200 watts and is changed from [watts] to [Btu/hr] with the conversion constant of 3.413. (10) The amount of heat transferred from the refrigerant to the airvapor mixture in the condenser is calculated as the sum of the heat transferred at the evaporator and the power added at the compressor. This assumes no heat losses in the refrigeration cycle. (11) This block of equations calculates the properties of the air-vapor mixture as it passes the evaporator, corresponding to equations (1) and (2) of the Basic Equations section. The amount of heat transfer in [Btu/lbm] is calculated using QOUT and the density of the air in region A. A value for the enthalpy of the mixture coming out of the evaporator is then found by subtracting this amount of heat transfer from the enthalpy of the mixture going into the evaporator. ENTHB = ENTHA - QEVAP [Btu/lbm] The mixture coming out of the evaporator is assumed to be at saturation. The PSYCKT subroutine is called to get the density of the mixture in region B and an average density and mass air flow in the evaporator is calculated: DENSEV = (DENSA + DENSB)/2 MDAEVA = CFM x DENSEV x 60 [lbm/hr] Using this average mass air flow, the heat transfer and resulting enthalpy is recalculated. To compensate for the amount of energy lost due to the condensed liquid which is drained away, a constant fraction of the enthalpy is subtracted. ENTHB = ENTHB -.0085 ENTHB The enthalpy of the condensed liquid is assumed to be the enthalpy of the saturated liquid at the same temperature. Since this value is very small relative to ENTHB this is an adequate approximation. The PSYCKT subroutine is called to calculate the remaining properties of the air-vapor mixture in region B. (12) This block of equations calculates the properties of the air-vapor mixture as it passes the condenser, corresponding to equations (3) and (4) of the Basic Equations section. The logic is identical to that of the evaporator 65

blcck except that here heat is added to the mixture and there is no condensed water to deal with. Also, the specific humidity of the mixture is assumed to remain constant over the condenser. With this value of specific humidity and the calculated value of enthalpy the PSYCKT subroutine is called to get the remaining properties of the air-vapor mixture in region C. The next section of the program constitutes the application of the four constraints discussed in the previous section of this report. As these constraints necessitate the change of the refrigerant temperatures or pressure, new values for the heat transfer quantities are calculated and used in blocks (11) and (12) above. (13) The control integer ISS remains zero until the high side pressure reaches the maximum value permitted, it is then set equal to 1. When this occurs, the only constraint checked for the remainder of the drying time is that the temperature of the air-vapor mixture going into the dryer is no higher than 5~ less than the condenser temperature. (14) If either the enthalpy of region A or the temperature of region B is higher than the stipulated design values, the properties of the refrigerant are changed. When the enthalpy in region A becomes higher than the designed value the evaporator temperature is increased by an amount proportional to the enthalpy difference. HDIFF = ENTHA - HMAX TADDE = 1.7 x HDIFF The constant of proportionality of 1.7 is determined from the McQuay data curves, samples of which are in the Appendix. When the temperature in region B becomes higher than the designed value, a check is made on the condenser temperature to ensure that it has correspondingly increased. When either or both of these changes in refrigerant temperatures is made, the heat transfer values are recalculated. The constraints and subsequent alterations of this section are checked three times before the program proceeds. (15) The temperature of the air-vapor mixture as it leaves the condenser can never be greater than 5~ less than the temperature of the refrigerant in the condenser. If an intermediate value violates this constraint the heat transfer in the condenser is reduced to the point where the constraint is met. QADD = QADD - QDIF 64

Since this reduction in heat transfer will increase the temperature of the refrigerant in the cycle, an incremental temperature difference chosen by the user is added to the temperature of the condenser and evaporator. CONDTP = CONDTP + CONADD EVTP = EVTP + EVADD Using these new temperatures the heat transfer in the evaporator is also recalculated. (16) The final constraint is to limit the high side pressure of the refrigeration cycle to a maximum value chosen by the user in conjunction with the capabilities of the compressor and refrigerant being used. When the pressure reaches this maximum, the control integer ISS is set equal to 1 and the heat transfer values and refrigerant properties that exist at that time remain constant for the remainder of the drying time. (17) Now that all the air-vapor mixture properties and refrigeration cycle properties have been made compatible with each other, these values are printed out in a temporary file which can be listed out when the execution of the program has terminated. These values are printed out every.3 min of real time until the maximum high side pressure has been reached. (18) The DRYER subroutine is called to perform calculations representing the drying of the clothes based on the inlet conditions and returns to the main program the enthalpy and specific humidity of the air-vapor mixture coming out of the dryer in addition to the new properties of the clothes load. (19) The PSYCKT subroutine is called to calculate the remaining properties of the air-vapor mixture. The mass flow of the dryer exhaust and drying rate of the clothes are also determined. (20) Every ten times the DRYER subroutine is called (or every 3 min), the properties of the clothes, the air-vapor mixture in the three regions, and the refrigeration cycle are printed out. (21) A maximum of two arrays are filled with the values of two properties of the system to be plotted against time. The properties to be plotted are set equal to FUNC1 and FUNC2. In the flow diagram these properties are arbitrarily chosen to be TDBA and CTEMP. A value is put into each array every time the values of the system are printed out, or every 3 min. 65

(22) The moisture content of the clothes is checked and the simulated drying system keeps operating until the water content of the clothes is less than or equal to 35 of the dry weight. (23) When the clothes are determined to be "dry" the final values of the system are printed out and the GRAPH subroutine is called to plot the graphs that the user has chosen. 7. DISCUSSION OF PSYCKT SUBROUTINE The subroutine performs psychrometric calculations to determine the properties of the air-vapor mixture. The properties used in the subroutine are: TDB = Dry bulb temperature [~F] TWB = Wet bulb temperature [~F] HUS = Specific humidity [Ibm H O/ lbm dry air] DENS = Density of the mixture [Ibm dry air/ ft ] ENTH = Enthalpy of the mixture [Btu/lbm dry air] Other properties used in the calculations but not input or output are: PW = Water vapor pressure at wet bulb temperature [in Hg] BARO = Total pressure of air-vapor mixture [in Hg] PVA = Partial pressure of water vapor in air [in Hg] The equations used in this subroutine were obtained from Richard Fanson of the Laundry Group of Whirlpool. These equations were used in the form in which they were obtained except for modification of the value used for the enthalpy of evaporation due to temperature variations. The subroutine includes four cases in which one or two of the properties are known and the others are calculated. These cases correspond to the value of the control integer NJT. 66

MAIN PROGRAM REAL MOAB/ MDACOA, //ODAC, AM'O/Trj MDAA MDADA, MfDEVA D/MEs/OA/,J X(3), XDoT),, /DR(25z) ~Y(25z), yz CZS),., | COM ON\/ X,001-j x A4J N9 AI j14PJ/CCO/\XC,, J TADB, /7T, T) 7-9 P e COA1HT-fF, =,CPA W.-rZ >W1 L, PZJ FZ CPC IG1I MMOIaT = C M o/ST x -OAO r C- 5 | WATT' = o TDA'= TD, TBSA = 7T,: | CY |,L PsyCA<T LSOB oUT7A/ ( 7TDSA, -rWBA1 USAD)N AT EwA,9 I | MP AA = CFM X DEASA > 0- (o | CAC- DRYE/ SC)S/ROUTrJe (O. JDBA,-TrSA,-rl/m tCT~PMMOQl 5TMDAA, XX, y, CQLOAD, WATTS, EAHC, gS C)., ~1,.... /,/=0. -.Z'" =O /P,.. = Cj /,,'Z. =0 Z'P7'= o i 0?0T 30 (3 00-P), 2-/.O(E.VTP-35-) + 9 | Powu) -.35 *. 06(v VSSp- 0 v)7r(g*P - 4.3-S-) ^ /030 + /.79 (ksP-(a) | %,

ouJ zoPO O =I ZOO y4o PO Li = woC),,. 3.4/3 9:ADD = QOOT -- PocAS () — EVAP = -QOUT/MDAA ~Al T- H = E/'- 7ZiA - ~G VA P CALL. P-y/KT SODeoO'T// (TTDeC,8)7-r0 T6, /5, D~&5A,/ ENT/B4) 4) DEAiSS/E\ (DECSA~ DECA2)/2 MDAEVA - CFM' DEAJEV Yx (0 EVAP- COUT/MDAEVA ETAr - EA/ rT/A -e EV4/ P EN'r/./ s E tr/'/ - oC, &ENr..E, C AL- PsycKT oEu? UTr/A E (TDE8, TEB, iu5 B, DENS6, ENA/ RB 4) o, MAP =A cFrM y 0 AM EAlC, Yo QCOA./I - Q AC C/ MOA % ENTT C'- E-TJ T -+ -QCOAIO HDe/r oC (C'v, + CALc PA\/L C T BS PORUT/~E (T-rDC,T)S C, HP 5C, DEC A/ EfCA/THC, 2) DE ACO = ( DetO5 R +- D E/SC)/ z M oACoa a- C FM v DE Sc/O Y FO QCr D - Q A D /MDA C.o/ EAjTbC _.- ^TIa + cQCoA/.O CAu- PsyCK/T s oB ouT/AIl (DBC, r76C/4CD A/C, CENT-C) Z) MDAC = C FM' DEo/SC (O> __ -- 0S = 0 T-MAX C/OADTP-O TR - 7C h 77 AX as, ~NOX s ~ ~~~~~~~~~~~~~~~~~~~~~J

|mAC = CFM K oeASC so I~-<^ ~EN-TA -L HMAX A/0~ TDB T3M-AX Yes TI I 77 T >_ — < TklA i kI MAX - HDIFF = ETV/7A - 14MAX Yes o rTADD = /. 7 X HOIF"F E'V TP = EV TP/* + T4A D OE VIES COQA/MI/A- CO2z/V/ F TA=OC L- < COhiOTP CoAlM/^/ = D Ies I sp= 5(CoDrP -'")/.25A (_

A A Y= CQoNrDP-5 L 7-DSc TMAX YE rD8Cc =7M AX CAL L Psc rT SoBeO 3 soT/ (71D0, 6W8C, W/SCC, OESSC, EASTAC 1,3) aIDAC1 = CFM ) O(E/SC Yw ( EIVD/F = ENA1h'C - eAlTVrCl QDIF= E^AD/O K (MDAC M DACl)/Z O QADD OA= D - QRD I COA/DT-P COAJD -P + CON^ADD | 5P (- (CO +A"T P - s )/. Z lo — < 5 P SHSPMA- xSP=PsPMAX yes i ___ E'VrP -EVTp EVA | | E A | qo-OTr = 30 (30Co -.SP) +.Zl 81, SOO (,(VTP -3 G) + 9Z 00 L-=L+1 ^-

L= 0 WRITE (ro7a, CQohtrP. -roe, eaEMNTA, EVrP, Qoour, QADO, TlME) ( y) —-\ CALL eYEeR Sd8eD(O-TtAe(.)rD6<i,7;T&<!,7T/M~,<CTEMPMMOtsr,MDAC7,T48A vsA,<4.OAD, WAv-T-SC, AIT-C4 ) CALL Pf:^rK- S V8gDVe-r1,(7-aA, T-rWBA, 1US4,,OEA1SA, -N-r,#A, MDAA = CFM x OENSA Y -0 DDeATE X ODOT( l) A/'~A/^J ~~L? I ~~I L /O I\E Wet-rE7 (77/ME,gMM0QrT;CTeMfiTQDc, clusg, T08A, T ysAT08a^ S sD,DtAre^AoDI QooT,ecoAD)A/orP Ev ZP PL =0- - = 1 - F NCO0NM = CTEMP PCok (rPr) =T/^t yz (rPr- )= CUIJC | ^^/<:; = TD6A A/PP?- = VPP? + AfIpT) =^/=/CI i ~~~~NPP1 = n/PPl+1l

-- MMCJOIST/CLOAD ^.03 ~^ —--- j AL PQyCKT' SOBoor/' TIG (eOArTSA,W.'5A, O'NSA ENT7A,) | yes. i ro

Case 1. The dry bulb and wet bulb temperatures are known. Since this case is used in the program with only ambient conditions, 1050 Btu/lbm is used for the enthalpy of evaporation which corresponds to 76~F. Case 2. The enthalpy and specific humidity of the air-vapor mixture are known. This case is used with the mixture as it enters the dryer. Therefore, the value used for enthalpy of evaporation is 1035 Btu/lbm which corresponds to 103~F. Case 3. The dry bulb temperature and specific humidity of the mixture are known. Since this case is used with the dryer exhaust the enthalpy of evaporation used is again 1035 Btu/lbm. Case 4. The enthalpy of the mixture and the fact that the mixture is at saturation are known. 8. DISCUSSION OF GRAPH SUBROUTINE This subroutine prepares the information necessary for calling on a subroutine stored in the MTS Library which prints out plotted graphs. Y1MAX and Y1MIN are the maximum and minimum values of the variable to be plotted on the y-axis of the first graph. Likewise, Y2MAX and Y2MIN are the maximum and minimum values for the second graph. If this program is run on a system other than that of The University of Michigan, IPLOT must be set equal to zero and no plots will be produced. 9. DISCUSSION OF DRYER SUBROUTINE The DRYER subroutine package which includes four subroutines was used to simulate the actual clothes dryer portion of the system. This package was previously developed by Whirlpool. The meanings of the symbols used as arguments of the subroutine are given in the main program discussion. The control parameter must be set equal to zero to initialize a new run and is thereafter set equal to 1. The DRYER package controls and increments the time parameter. One alteration of the subroutine was necessary to use it in a ventless dryer system. An additional call to the package's own psychrometric subroutine was needed to calculate new properties of the air-vapor mixture entering the dryer every time the DRYER subroutine is called. In a conventional dryer the inlet air properties remain constant but in a closed-loop system these properties constantly change. 73

PSYCHROMETRIC SUBROUTZIE S o8 eoT/E SyC Kt (rDB, TS, VS, DeS/, ENT/, 7-r) | I4eo =~.2 9.' | =./6&943 x /DZ 7r -2 X =(Ex.8-,.//5,099/1 X/0-4(T-r2,0/00/oo/ 3 TS7 -. 7/40o957) | Pc = Exp'[(-a-x)/-. /4/739 1 PVA p wpJ - (A Ro - YPXrB -Tu)S) /(Z 80 - 3a 7W|S) <-, | = hS =, Z4Z PI/A(BAReo-PVA) Cz~,PVA~8Aot, |EH/r -.7.4 rDT3 +,#As /oso c.^ rDB) | | T7S = (El/r - /o3S- /us)/(.24 +,44 A406) ( d' —- PVA =.u., BA o/(. (, 2. +.4 a u4)l CL co 2900 Y PVA + TDB x 13tARo | CZ=:.3 P/VA +,SAA.o C3 2700 + rD. Tu)B 0 d —@ c

-- eXi = -T2 a | S= 2o /6z 43 x/0o-z T''S -1 X - = [8-.//5o 98/ x/0o'' Tt'B. +,O/OO/ 283 3 Tr -. 7/140956 7?' Pw = EX P [(6 }X)/. /4/ 73497 FX=2 PWU- (CZ- CZ7-Ja)/(c3.- -.3 Yr 8) FX1= PW)/. /41/7349 [(.,D(493~XD /5 T6 t3 -.//.o910/ x/o-4u TrW 70.oo84/5 x/o-L) / X +. /L943 /ot ] - (C3 V. C 2-.3 C?C) /(C3-3 z. a)|T~JI - x - FX/FXI I,T0 \J1 1- V |-' - ~t3a ~.o5 |O5~ ( - EAT^ e." Al TOe+^US(/03^ u.44 7- --- t-^ 7-oe em -rer/(oo 00 5?47<O G6 EA —-. Z 90 o /i)| 5 *S =. 00/40973 9 EX P (C o3 4 0 (S7 s TD7 ) | |r B 7-T0 | R2A -B A RO ~ SVS/.La;} / ^) I I @ ( ) —"I ODENS (GA Ro,-PVA)/[.754/(460o+TD8))j T,I~ ~.- ---- ~,, ETO

GRAPH SUBROUTINE SU (3o orAE GRA PH (/rPLoT /PPI, /PPJ e, y j, y 2) P/^MEA/oSoJ so//oR (25), y (z2), yz2('S), z/rMACE(oo) IPL7 o > 0 CALL PLO-rT14 -ATS LIB/SAR SOQBeoUTr/4E.ALPCLT >1 No )/ eOEr8/sQ7 — 4,','__,'"',... —— ", _ ZMAX =/5", yZM/I /= ao CALLe PLoTi4 -/TS L46AR/ SVoUT',e/,VE'

10. EXPERIMENTAL RESULTS To evaluate the possible capabilities of a clothes dryer system of this type and to obtain some initial design parameters the effect of several parameters on the operation of the system was investigated. The input data which were kept constant throughout these runs were: Ambient dry bulb temperature TDB = 72~F Ambient wet bulb temperature TWB = 60~F Initial clothes load temperature CTEMP = 72~F Initial high side pressure of refrigeration cycle HSP = 290 psig Initial evaporator temperature EVTP = 45~F Maximum high side pressure HSPMAX = 400 psig Designed air temperature over the condenser TBMAX = 60~F Temperature added to condenser in constraints CONADD = 3.0~F Temperature added to evaporator in constraints EVADD = 0.38~F (1) In Figure 5 the effect of moisture content and load size on drying time is shown. The remaining input data were: Volumetric air flow CFM = 220 CFM Designed air enthalpy over the evaporator HMAX = 40 Btu/lbm while varying the moisture content: Dry weight of clothes load CLOAD = 9 lb 77

O4100 Iz w F- e12 0 - 0 O C-) w 0m. 075 9 cn 6 ~ w I-, / — CLOTHES LOAD z j --- MOISTURE CONTENT 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 DRYING TIME (HRS) Figure 5. Clothes load and moisture content vs. drying time. 78

while varying the clothes load: Initial moisture content of clothes CMOIST = 85% (2) The effect of the volumetric air flow in the closed-loop on the drying time required is shown in Figure 6. The remaining input data were: Designed air enthalpy over the evaporator HMAX = 40 Btu/lbm Dry weight of clothes load CLOAD = 9 lb Initial moisture content of clothes CMOIST = 85% (5) The effect of increasing the designed air enthalpy over the evaporator is an increase in the amount of time into the drying period before the constraint which increases the evaporator temperature begins to operate. Since an increase in the evaporator temperature causes an increase in the capacities of the system, the longer the time before this increase begins the longer the drying time required. This effect is shown in Figure 7 where the remaining input data were: Volumetric air flow CFM = 220 CFM Dry weight of clothes load CLOAD = 9 lb Initial moisture content of clothes CMOIST = 85% The following is an example of the computer simulation output. The input data are given in the printed output. The results given here are those used to compare the simulated system to the other systems in Figure 2. Following the program output is a listing of a temporary file containing the values of properties in the system before the maximum high side pressure was reached. It is evident that with the refrigerant, compressor, and resulting maximum pressure used here, the fluctuations in heat transfer at the evaporator and condenser occur during only the first few minutes of the drying time. The capacity of the condenser during this time is shown in Figure 8. The capacity drops as the condenser temperature increases and then increases as the evaporator temperature starts increasing until the maximum pressure is reached and the capacity remains constant thereafter. The capacity of the evaporator follows a similar pattern. 79

300 g 250 u200 0.7 0.8 0.9 1.0 DRYING TIME (HRS) Figure 6. Air flow vs. drying time. 80

45 I / 40113 w 35 0.90 0.95 1.00 DRYING TIME (HRS) Figure 7. Enthalpy vs. drying time. 81

16,000 #sof-o z o 15,000 LQL z 0 0 14,000 0.01.02.03.04.05.06.07 TIME (HR) Figure 8. Condenser capacity vs. time. 82

**m** DRYER SIMULATION ***** DRYER INPUT DATA AMb. DRY d-ULB TEMP. = 72.0 F AMB. nET bULB TEMP. = 60.0 F AIR I-LOw = 220.0 CFM INITIAL CL01HES TEMP. = 72.0 F INITIAL CLO1HES LOAD = 9.0 LBS INITIAL MOISTURE COJTENT = 85.0' REFRIGERATION CYCLE INPUT DATA CONDENSER EVAPORATOR INITIAL PRESSURE = 290.0 PSI INITIAL TEMP. = 45.0 F MAXIMUM PRESSURE = 400.0 PSI DESIGNED AIR TEMP. = 60.0 F DESIGNEO AIR ENTHALPY = 40.0 BTU/L3M DELTA TEMP. = 3.0 F DELTA TEMP. = 0.4 F THE RESULTS FOLLOW TIME MMOIST CTEMP TDBC HUSC TUBA HUSA TD[0 HJSB DRATE QAUD QOUT CONDTP EVTP o 30.05 7.288 82.5 141.0 0.0142 87.0 0.0238 o7.9 0.0142 -8.1 15074.2 11240.0 148.5 54.0 3.10 O.85 2 92.5 150.0 0.0221 91.5 0.0333 80.8 0.0221 -9.2 15758.1 12228.7 155.0 61.8 0.15 6.400 96.0 150.0 C.0261 10.17 0.0314 85.7 0.021l -9.1 15758.1 1222. 7 155.0 61.8 0.20 5.948 97.2 150.0 0.0278 101.9 0.0390 87.5 0.0278 -9.1 151 155. 12228.7 155.0 61.8 0.25 5.496 97.7 150.0 C.0285 102.3 0.0396 88.2 0.0285 -9.0 15758.1 12228.7 155.0 61.8 3.30 5.045 97.9 150.0 C.0288 102.5 0.0399 88.5 0.0288 -9.0 15758.1 12228. 7 155.0 61.8 0.35 4.593 98.0 150.0 C.0288 102.6 0.0400 88.6. 0288 -9.0 157 58.1 12228.7 155.0 61.8 0.40 4.142 98.0 150.0 C.0289 102.6 0.u400 88.6 0.0289 -9.0 15758. 1 12228.7 155.0 61.8 0.45 3.691 98.1 150.0 0.0289 102.6 0.0401 88.7 0.0289 -9.0 15758.1 12228.7 155.0 61.8 0.50 3.240 98.1 150.0 C.0289 102.6 0.0401 88.7 0.0289 -9.0 15758.1 12Z28.7 155.0 61.8 0.55 2.789 98.1 150.0 C.0289 102.7 0.0401 88.7 0.0289 -9.0 15758.1 12228. 1 155.0 61.8 0.60 2.338 98.1 150.0 C.0289 102.7 0.0401 88.7 0.0289 -9.0 15758.1 12228.7 155.0 61.8 0.65 1.88 7 98.1 150.0 0.0289 102.7 0.0401 88.1 O.289 -9.0 15758.1 12226.7 155.0 61.8 0.70 1.454 100.5 150.0 C.0275 104.9 0.0377 87.2 0.0275 -8.2 15758.1 12228.7 155.0 61.8 0.75 1.069 106.2 150.0 0.0263 1)10.1 O.O050 85.9 0.0263 -7.1 15758.1 1222d.7 155.0 61.8 0.80 0.752 114.2 150.0 0.0255 117.4 0.0324 85.0 0.0255 -5.6 15758.1 12228.7 155.0 61.8 0.85 0.513 123.4 150.0 C.0250 125. 0.0300 84.4 0.0250 -4.0 15758.1 12228.7 155.0 61.8 0.90 0.348 132.0 150.0 C.0250 133.6 0.0283 84.5 0.0250 -2.6 15758.1 i2228.7 155.0 61.8 0.93 0.269 136.8 150.0 0.0253 138.0 0.0277 84.8 0.0253 -1.9 15758.1 12228.7 155.0 61.8

SLIST -A I TDIBC CONDTP TOBB ENTHA EVTP WOUT QADD TIME 2 103.1 127.5 38.6 26o.4 45.0 11748.0 15474.1 0.0 3 120.9 127.5 51.5 33.7 45.0 11748.0 15474.1 0.005 4 123.1 130,5 54.1 35.1 45.4 11473.4 15182.6 0.010 5 126.1 133.5 56.3 36.3 45. 8 11198.8 15144.1 0.015 6 129.1 136.5 58. 7 37.7 46. 1 10924.3 15091.8 0.020 7 131.0 136.5 60.8 39.1 46.1 10924.3 14908.7 0.025 8 132.1 139.5 63.1 40.6 46,4 10616.6 14470.8 0.030 9 135.1 142.5 64.7 42.0 48.8 10806.9 14665.6 0.035 10 138.0 145.5 66.3 43.5 51.3 11014.8 14865.5 0.040 11 141.0 148.5 67.9 45.0 54.0 11240.0 15074.2 0.045 12 144.0 151.5 69.4 46.6 56.6 11482.0'15292.0 0.050 13 147.0 154.5 70.9 48.2 59.4 11740.1 15519.8 0.055 END OF FILE co

APPENDIX (1) Figure A-l is a schmatic diagram of a system where the refrigeration system is used in an open-loop configuration. Such a system operates more effectively than a closed-loop system because the air entering the condenser is ambient air at approximately 50to relative humidity rather than saturated air. The air into the dryer is also at a "dryer" state and more able to absorb the moisture of the clothes. An additional feature of such a system would be the possibility of using either the cool air from the evaporator or the warm humid air from the dryer to cool or humidify the room. (2) Calculation of the energy lost in draining the condensed water. This calculation is made using the data from sample output of the simulation program. General equation: h + Qm hAV + B + ( - )h h A A B A B Mass flow of air: m = V/v x 60 = 100ft3mi x60 = 810 bm/hr 14.8 ft3/lbm x 60 Heat transfer: 12228 Btu/hr m = 810 lbm/hr = 15.1 Btu/lbm ENTHA = 69 Btu/lbm, HUSA =.0400 ibm H 0/lbm air Ignoring the drained water we get: ENTHB = 69 Btu/lbm - 15.1 Btu/lbm = 54 Btu/lbm 85

TO OUTSIDE AIR VALVE DRYER TO ROOM EXPANSION VALVE _- _X ROOM AIR EVAPORATOR CONDENSER DRAIN COMPRESSOR Figure A-i. Schematic diagram of open-loop system. 86

Since this is at saturation: HUSB =.0295 lbm H 0/lbm air Assume the enthalpy of the condensed liquid is the enthalpy of the saturated liquid at 880F. An approximate value for the energy of the liquid is: ( A - C )h = (.0400 -.0295) 56 = 1.1480 Btu/lbm A B I In the computer program, this is assumed to be a constant fraction of ENTHB: 1.1480/54 =.0212 The simplified approximation is then: ENTHB = ENTHB -.0212 ENTHB The value in the program listing is lower than this since it was calculated using lower temperatures. This higher value will be useful in lengthening the energy build-up time. (3) Calculation of the amount of air bleed at the dryer exhaust necessary for the system to remain at a steady state operating condition after the maximum high pressure has been reached. The values used are taken from the sample simulation output. Dryer exhaust conditions: TDBA = 102.6~F, HUSA =.040 Ibm H20/lbm air ENTHA = 69 Btu/lbm, v = 15 ft3/lbm relative humidity = 87o 87

Ambient air conditions: TDB = 72~F TWB = 60~F ENTH = 26.5 Btu/lbm, v = 13.6 ft3/lbm relative humidity = 500o Amount of energy that must be dissipated: QADD - QOUT = 15758 - 12228 = 3,530 Btu/hr Necessary air bleed: 3,530 Btu/hr bleed =(69-26.5)Btu/lbm = 8 lbm/hr Vbleed = 83 lbm/m x1 l 20.8 CFM bleed / 60 Approximate amount of moisture this would add to the room: (HUSA - HUS) m =.032 Ibm H20/lbm air x 83 Ibm air/m = 2.65 lbm H O/m. (4) Figure A-2 is a copy of the compressor curves used in the simulation. Another compressor could be used in the simulation as long as curves such as these were available from which to derive the necessary equations. (5) Figure A-3 shows the McQuay data curves relating evaporator temperature to the enthalpy of the air over the evaporator. (6) Figure A-4 shows the McQuay data curves relating condenser temperature to air temperature over the condenser. (7) The following is a compiled listing of the main program and subroutines PSYCKT and GRAPH. A sample of a computer output plot and data file are also included. 88

Revised i..-4- 6 Index io. -.:..3 --— tr^-itrj. TEiCM H PRnODUCT" iS CO. "EEE L WMoa2 o.:- _ +-: —::e_ TE CU.s-, M,4,c(.:-;,:~ E1 16. 13 il::::Xj^::=lli I1I!:1 qI3rll:qllliljii ^F::qz~Fllrir^^^^^:/tS~ _~~F~~eU DATA SHEET 3.t tS-i: - ei i ^l 3 4S ^ ^S 1 i i I ""' -II I — N5. m3';0~~-~~~~-ttdi~~~~~~..:;.': - i r. t n' i j!-.-'. —S - -^ - ff — r i A..f' - J " "., - j-, —t _ - ~!', s I' i " l _,' i L' t!-; mR=F. b A-. oeBo e 3 |~ i t' |! I [I l i I' ^ -j ^ S j: ^-r 1-. _ 25 30 33 40 45 50 __. _-_4-i- 4'. r: ^ i n i I t:1 ^ - * * ~' t r'p j i t p''T Tti^ -r'' i!: N " ^ 7 g'f ~ x y~ m T T ^ I~ t ^ ^' j!: ^3 1:'CITCi_ Jfcr89~

*SGn~aTIO aOosaaduloo ~-'v aanDIAT 4~~i1i~-ci? t~ r i- t~- I +-'rtlt-~vt- 7 t t" 7 —- +' 1T T -1 i4.- i 1 1 ~' ii j t i r. iL1 ~;-lj -t i tt ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ C' t44 i 14 -i Lif;; i ~c-? -!-r ~~~'t-t-t~~t-~~ t~ir ic: i -t-ir I I 4 i i4- t ri::7 i~~~~t:;:::~~~~~.:; 1' F F, ~ ~ ~ ~- ltT:17 F 1 tt ~(-"-ti-;- i c; c4 j-t 4 -I 4 t tr A~~ ~ ~ t I I F t; L ~ 1 Cc 7 4~i!,?;~; ri i~~~~~Ert-:r 44A+~~~~~~~~~~~~~~~~~ 4- L'4 4-t~~~~~~~~~if-i -i~ -t-t eceit~i He~L~ t(: ~ ~c-it: ~~;~~~~~~~~~~~~~~~~~~~i t f ~ ~ ~:i " i-e 1i I I L I Li ~ ~ ~ ~ li i t ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ t ~, tt r~ii l~r~. c~-~IFTKS~i~4FJ ~ ~ r j f7 4, t 4; t-~~tc~~~ f-~~~ ~-t i A T~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i _ bl~~~~~Cr~~~~J I;:' ( ~-Y~~~~~~~~~~-(-icC~~~~~e; U Ckirreci~~~~~~~~~~e~~c~~e ~~~~e.~! ii'-C"~~~~~~~4 +f ~t? 144 4-~;t

*S;9AJTIn TossaadmiUoo *+[-V 9-nSEj: r-t ~~~~..~~~~~~~~~4 ^ wil {tp ^ -:piwrpw-^^ a~w^ SS S ^ SE S ^^^^^^^^TO^ ^^r^ ^ ^1^^ ^ ^^^S^^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~. ~ o<^^^^^^^a^^40^^4 —^^ —^^^IS| ^^ ~g i ~iP fe^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~7 ^f^ — l -- -^i-^ -|^^ ^^^^ ___^ ^^^^H^^_ _^_^ni^^-*7-^-^- ns^^^^-^^-^^^iu-^I^^M^^M^E^II ^^ill~~~~~~~~~~~~~~~~~~~~~~~~~~~~~aftt^^~~~ s-^ ^^^^^^ly^^^^ll^^^- I ^|,|^^^.~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~itt~~~~tt~~~~~t^^^~i??eT^^^^-~-^-^^^^^^^ S~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i^SM~~~~~~~iiI^~~Bul i~~lif^^ ^-7Nj...... 44 6f, -^^ ^^^^ ^-^^^ —^^^^^^Xll!^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~4-T ^f^-^ ^^^fe^^^^a^-X-^il~~~~~~~~~~~~~~~~~~~aM^~a ^^^ ^M~~~~~~~~~~~~~~~~~~~iff-a^^^i-^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ — Ill iri?^ -7^^4 —^^^^^-^IL~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~a^-~~4 1?7 4 i N, t4TT 1 tti~~~ijN?^^ft ^^^^^:^ir^^^Rio Ja i+jiS l *^^ —^^7^ —-~ —^^ —— ^M —_.-__l^^M^:t^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *,,^^ ^-l,,^i^'/U.^^;:: ^: i iffi ^ i.,,,11.;!';;- wb'-i:^ ^ii~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~t;::;:::;;;: ^-^ ^ *sF^^^^ ^^ ^'6~^'"'"'"""~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ^ ^ W^H?^ ^ ^ —-TT-~^ ^T-^ ^ ^ ~ —ta^^^^~~~~~~~~~~~~41-- -— fci^_.A. ^^*-'^ **:;:^; ^::::: ^;i;;.~ *;;* * *;'*'- -- -- -- -- — I....... — ** -'*; *"*"'"y ^^ -^;s^; *:;.;.*.5k'r*;.*I- *s.?'"^ ^ —-* —-^-;-rS- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ B,3l||ii-:||^j|J lli~~~~~~~~~~~~~~~~~~~~liiiM~~~~~~~ilJ~~~~~iS^^I^^^l~~~~~~~~~~t^^^^^^tt^ tt~~~~~~~~~~~~~l s^^^^^~~~taeTi ^'t ^^'^ "I'~+-*Ti~? s Ty -^ -4~ -^.: ~ -0: ~ -i~S s S;ii S -ii -^ l I -rl^-l~tri 1,; ^!pp4^1^^4ti1 }a(^ J L,~i ^ lt ~ ^ ^ ]^ O 4|j^ l~ i^ g TI ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~?'^^

FORTRAN IV G COMPILER MAIN 05-10-71 16:28.00 PAGE 0001 C************************ DRYER SIMULATION ************************ C** ** C** THIS IS A SIMULATION OF A CLOSED-LOOP CLOTHES DRYER WHICH ** C** PASSES THE DRYER EXHAUST AIR OVER THE EVAPORATOR AND ** C** CONDENSER COILS OF A REFRIGERATION CYCLE IN THE FEEDBACK ** C** LOOP. ** C**** C** THE PURPOSE OF THIS SIMULATION IS TO AID IN THE DESIGN OF ** C** SUCH A CLOTHES DRYER SYSTEM. ** C**** C** THE FOLLOWING IS A LIST OF THE INPUT AND OUTPUT PARAMETERS:** C** INPUT: ** C** TDB - AMBIENT DRY BULB TEMP. (F) ** C** TWB - AMBIENT WET BULB TEMP. (F) ** C** CLOAD - DRY WEIGHT OF CLOTHES LOAD (LBS) ** C** CMOIST - MOISTURE CONTENT AS PERCENT OF DRY ** C** WEIGHT M(% ** C** CTEMP - INITIAL CLOTHES LOAD TEMP. (F) ** C** CFM - VOLUMETRIC AIR FLOW (CFM) ** C** HSP - INITIAL HIGH SIDE PRESSURE (PSI) ** C** EVTP - INITIAL EVAPORATOR TEMP. IF) ** C** HSPMAX - MAXIMUM HIGH SIDE PRESSURE (PSI) ** C** HMAX - DESIGNED AIR ENTHALPY OVER EVAPORATOR ** C** ( BTU/LBMI ** C** TBMAX - DESIGNED AIR TEMP. OVER CONDENSER (F)** C** CUNAUD - TEMP. ADDED TO CONDENSER WHEN AIR ** C** TEMP. GETS TOO HIGH ** C** EVADD - TEMP. ADDED TO EVAPORATOR WHEN CONADD** C** IS ADDED TO CONDENSER ** C** IPLOT - CONTROL INTEGER PLOTTING: ** C** 0- NO PLOTS ** C**I - ONE PLOT ** C**2 - TWO PLOTS ** C** ** C** OUTPUT: ** C** TIME - TIME IN HRS ** C** MMOIST - MOISTURE CONTENT OF CLOTHES (LBS) ** C** CTEMP - CLOTHES TEMP. IF) ** C** TDBA,TDBBTDBC - DRY BULB TEMP. OF AIR ** C** LEAVING THE DRYER, LEAVING THE ** C** EVAPORATOR, AND ENTERING THE DRYER, ** C** RESPECTIVELY lFa ** C** HUSA,HUSB,HUSC - SPECIFIC HUMIDITY OF THE ** C** AIR LEAVING THE DRYER, LEAVING THE ** C** EVAPORATOR, AND ENTERING THE DRYER, ** C** RESPECTIVELY (LBS WATER/LBS AIR) ** C** DRATE - DRYING RATE (LBS/HR) ** C** QADD - CONDENSER CAPACITY (BTU/HR) ** C** QOUT - EVAPORATOR CAPACITY (BTU/HR) ** C** CONDTP - REFRIGERANT TEMP. AT CONDENSER ** C** EVTP - REFRIGERANT TEMP. AT EVAPORATOR ** C**** C C 92

FORTRAN IV G COMPILER MAIN 05-1.0-71 16:28.00 PAGE 0002 0001 REAL MDAB,MDACON,MDAC,MMOIST,MDAA,MDAD,MDEVA 0002 DIMENSION X(3),XDOT(.3),HOR(25),YI(25),Y2(25) 0003 COMMON X,XDOTHAH,HHCAP,HCONV,XICRtTIUB, IT,THTVP,CONHT,EFF,CPAW,T2 1, WTCL,H2,E2,CPC,G1, QIN,T3DB 0004 NAMELIST/OATA/TDB,TWi,CFMtCTEMPCLOADCMOI ST,HSP,EVTPHSPMAX, I IPLOT,HMAX,TBMAX,CONADD,EVADD 0005 24 FORMAT('-THE RESULTS FOLLOW' ) 0006 25 FORMAT('-TIME MMOIST CTEMP TOBC HUSC TOBA HUSA TD 188 HUSB ORATE QADD QOUT CONDTP EVTP ) 0007 30 FORMAT(F5.2,F8.3,2F8.lt,F8.4,F8.1,F8.4,F7.1,F8.4,Fd.1,4F10.1) 0008 3 FORMAT ('-',15X,'DRYER INPUT DATAt) 0009 4 FORMAT('O','AMB. DRY BULB TEMP. = *,F5.19, F',5X, l'AMB. WET BULB TEMP. =',F5.1,' F'/IX,'AIR FLOW = It 2F5.1,' CFM',I4X,'INITIAL CLOTHES TEMP. =',F5.1,' F'/IX, 3'INITIAL CLOTHES LOAD =',F5.l,1 LBS',2X, 4'INITIAL MOISTURE CONTENT =',F5.1,1 %') 0010 80 READ (5,DATA) 0011 WRITE (6,8) 0012 8 FORMAT (' 1',20X,'***** DRYER SIMULATION *****') 0013 WRITE (6,3) 0014 WRITE (6,4)TDB,TWB,CFM,CTEMP,CLOAD,CMOIST 0015 5 FORMAT('-',15X,'REFRIGERATION CYCLE INPUT DATA') 0016 WRITE (6,5) 0017 WRITE (6,6) 0018 6 FORMAT ('O','CONDENSER',25X,'EVAPORATOR') 0019 WRITE (6,7) HSP,EVTP,HSPMAX, TBMAX,HMAX,CONADDEVADD 0020 7 FORMAT('O','INITIAL PRESSURE =',F5.1,' PSI',6X, 1' INITIAL TEMP. =',F5.1,' F'/IX,'MAXIMUM PRESSURE =', 2F5.1,' PSI'/IX,'DESIGNED AIR TEMP. =',F5.1,' F',6X, 3'DESIGNED AIR ENTHALPY =',F5.1,' BTU/LBM'/IX, 4'1DELTA TEMP. =',F5.1,' F',13X,'DELTA TEMP. ='tF5.1, 5' F') 0021 MMOI ST=CMOIST*CLOAD/100. 0022 WRITE (6,24) 0023 WRITE (6,25) 0024 9 FORMAT (5X,'QADD QOUT CUNDTP EVTP') 0025 WATTS=O. 0026 EVTPIN=EVTP 0027 TDBA=TDB 0028' TWBA=TWf 0029 WRITE (7,2000) 0030 2000 FORMAT (4X,'TDBC CONDTP TUBB ENTHA EVTP QOUT', 14X, QADD TIME') 0031 WRITE (8,2000) C INITIALIZE DRYER SUBROUTINE 0032 CALL PSYCKT(TDBA,TWBA,HUSA,DENSA,ENTHA,1) 0033. MDAA=CFM*DENSA*60. 0034 N=0 -0035 CALL DRYER(O,TDBA,TWBA,T IME, TEMP,MMOI STMDAA,XX,YY,CLOAD,WAITS EN 1THCHUSC) 0036 J=0 0037 L=0 0038 ISS=0 0039 NPP1 =0 0040 NPP2=0 93

FORTRAN IV G COMPILER MAIN 05-10-71 16:28.00 PAGE 0003 0041 IPT=O 0042 CONDTP=. 25*HSP+55. 0043 CONIN=CONDTP 0044 500 BTUEV=30.*(300.-HSP)+.281*800.*(EVTP-35.1+9200. 0045 POWW=(4.35+.06*(HSP-260.) )*lEVTP-43.5)+1030.+1.75*(H-SP-260.) 0046 IF (POWW-1200.) 520,520,510 0047 510 POWW=1200. 0048 520 POWB=POWW*3. 413 0049 BTUCON=BTUEV+POWB 0050 QOUT=BTUEV 0051 QADD=BTUCON C AIR IS NOW ENTERING EVAPORATOR. 0052 600 QEVAP=QOUT/MDAA 0053 ENTHB=ENTHA-QEVAP 0054 CALL PSYCKT(TOBBiTWBB,HUSB,DENSB,ENTHB,4) 0055 DENSEV= ( DENSA+DENSB) /2. 0056 MDAEVA=CFM*DOENSE V*60. 0057 QEVAP=QOUT/MOAEVA 0058 ENTHB=ENTHA-QEVAP 0059 ENTHB=ENTH8-.0085*ENTHB 0060 CALL PSYCKT(TD8f, TWBB,HUSB,DENSB,ENTHB,4) 0061 MDAB=CFM*DENSB*60. C AIR IS NOW ENTERING CONDENSER. 0062 QCOND=QADD/MDAB 0063 ENTHC=ENTHB+QCOND 0064 HUSC=HUSB 0065 CALL PSYCKT(TOBC,TWBC.HUSC,DENSCENTHC~2) 0066 DENSCO=(O DENSB+DENSC)/2. 0067 MOACON=CFM*DENSCO*60. 0068 QCOND=QADDO/MOACON 0069 ENTHC=ENTHB+QCOND 0070 CALL PSYCKT(TOBC,TWBCHUSC,DENSC,ENTHC,2) 0071 MOAC=CFM*DENSC*60. 0072 IF (ISS-L) 602,601,601 0073 601 TMAX=CONDTP-5. 0074 IF (TDBC-TMAX) 660,660,900 0075 900 TDBC=TMAX 0076 CALL PSYCKT (TDBC,TWBC,HUSCtDENSC,ENTHC,3) 0077 MDAC=CFM* DENSC*60. 0078 GO TO 660 0079 602 IF (L-1) 603,640,658 0080 603 IF (ENTHA.LE.HMAX.AND. TDB8.LE.T8MAX) GO TO 640 0081 J=J+1 0082 IF (J-3) 605,605,604 0083 604 J=O 0084 GO TO 640 0085 605 IF (ENTHA.LE.HMAX) GO TO 606 0086 HDIFF=ENTHA-HMAX 0087 TADDE=1.7*HDIFF 0088 EVTP=EVTPIN+TADDE 0089 IF (TOBB.LE.TBMAX) GO TO 620 0090 606 TAODDC=TDBB-TBMAX 0091 CONMIN=CONIN+TAODC 0092 IF (CONDTP-CCNMIN) 610,620,620 0093 610 CONDTP=CONMIN 94

FORTRAN IV G COMPILER MAIN 05-10-71 16:28.00 PAGE 0004 0094 HSP=(CONDTP-55.)/.25 0095 620 GO TO 500 0096 640 TMAX=CONDOTP-5. 0097 WRITE (8,6421 TDBC,CUNDTP,TD88,ENTHAEVTPQOUTQADDO,TIME 0098 642 FORMAT (7F8.1.F8.3) 0099 IF (TOBC-TMAX) 658,658,650 0100 650 TODBC=TMAX 0101 CALL PSYCKT(TOBC,TWBC,HUSC,DENSCtENTHCl,3) 0102 MDAC I=CFM*DENSC*60. 0103 ENDI F=ENTHC-ENTHC 1 0104 QOIF=ENDIF*(MDACl+MDAC)/2. 0105 QADD=QADO-QD IF 0106 CONDTP:CONDTP+CONADD 0107 HSP= (CONOTP-55.)/.25 0108 IF (HSP-HSPMAX) 654,654,652 0109 652 HSP=HSPMAX 0110 ISS=1 0111 CONDTP=.25*HSP+55. 0112 GO TO 655 0113 654 EVTP=EVTP+EVADD 0114 655 QOUT=30.* ( 300.-HSP) +. 281*800.*( EVTP-35.) +9200. 0115 L=L+1 0116 GO TO 600 0117 658 L=0O 0118 WRITE(7,642) TDBC,CONOTP,TDBB,ENTHA,EVTP,QOUT,QADD,TIME C AIR IS NOW ENTERING DRYER. 0119 660 CALL DRYER(1,TDOC,TwBCtTIMECTEMPtMMOISTtMDACTODBAHUSAtCLOADWATT I S,ENTHC HUSC) 0120 CALL PSYCKT(TDBA,TWBA,HUSA,UENSA,ENTHA,3) 0121 MDAA=CFM*DENSA*60. 0122 DRATE=XDOT(1) } 0123 N=N+ 1 0124 IF (N-10) 680,670,670 0125 670 WRITE(6,30)TIME,MMOIST,CTEMP,TTDBCHUSCTDBAHUSA,TOBB,HUSB, I DRATE,QADD,QGUT, CONDTP,EVTP 0126 IPT=IPT+l 0127 IF ( IPLOT-11 676,674,672 1128 672 CONTINUE J129 FUNC2=CTEMP 0130 Y2(I PT)=FUNC2 0131 NPP2=NPP2+1 0132 674 HOR(IPT)=TIME 0133 FUNCI=TDBA 0134 YI ( I PT)=FUNC 1 0135 NPP.=NPP 1+1 0136 676 CONTINUE 0137 1000 FORMAT(4F10.1) 0138 N=O 0139 680 IF (MMOIST/CLOAD-.03) 700,700,690 0140 690 CALL PSYCKT(TDBA,TWBA,HUSA,DENSA,ENTHA,3) 0141 GO TO 600 0142 700 WRITE(6,30)TIME,MMOIST,CTEMP,TDBC,HUSC,TDBAHUSA, TDBBHUSB, 1 DRATE,QADD,QOUT,CONOTP,EVTP 0143 CALL GRAPH(IPLOT,NPP1,NPP2,HOR,Yl,Y2) 0144 GO TO 80 95

FORTRAN IV G COMPILER MAIN 05 -10-71 16:28.00 PAGE 0005 0145 END TOTAL MEMORY REQUIREMENTS 001330 BYTES 96

FORTRAN IV G COMPILER GRAPH 05-10-71 16:28.06 PAGE 0001 0001 SUBROUTINE GRAPH(IPLOT,NPPlNPP2,HOR,YL,Y2) 0002 DIMENSION HOR(25),Y(25),Y2(25),IMAGE( 1500) 0003 WRITE (6,5) 0004 5 FORMAT (*'1) 0005 IF (IPLOT-1) 30,10,10 0006 10 YlMAX=150. 0007 YIMIN=80. 0008 CALL PLOT14(Ot6g6,I,9,IMAGEt1.,0.,YlMAXtYIMINt*',tHORtYltNPPIl,4,O 10) 0009 WRITE (6,5) 0010 IF (IPLOT-1) 30t30,20 0011 20 Y2MAX=150. 0012 Y2MIN=80. 0013 CALL PLOT 14(0,6, 66,9, IMAGE, 1.0.,Y2MAXY2MIN,**',HOR,Y2,NPP2,4,0, 10) 0014 30 CONTINUE 0015 RETURN 00 16 END TOTAL MEMORY REQUIREMENTS 001A76 BYTES 97

FORTRAN IV G COMPILER PSYCKT 05-10-71 16:28.07 PAGE 0001 0001 SUBROUTINE PSYCKT(TB, TWB,HUS,DENS,ENTH,NJT) 0002 BARO=29.5 0003 IF(NJT-2) 100,200,10 0004 10 IF(NJT-4) 300,400,1000 0005 100 B=.1626943E-2*TWB-1. 0006 X= (B*B-.12150981E-04*TWB*TWB+.010016283*TWB-.71409557)**.5 0007 PW=EXP((-B-X)/(-.1417349)) 0008 PVA=PW- B ARO-PW)*(T TB-TWB)/(2800.-1.3*TWB) 0009 HUS=.622*PVA/(BARO-PVA) 0010 ENTH=.24*TDB+HUS*( 1050.+.44*TDB) 0011 GO TO 600 0012 200 TDB=(ENTH-1035.*HUS)/(.24+.44*HUS) 0013 210 PVA=HUS*BARO/(.622+HUS ) 0014 CI=2800.*PVA+TDB*BARO 0015 C2=1.3*PVA+BARO 0016 C3=2800.+TDB 0017 TWB-O 0018 4 X1=TWB 0019 8=. 1626943E-2*TW8-1. 0020 80 X=(B*B-. 12150981E-4*TWB*TW8+*.010016283*TWB-.71409557)**.5 0021 PW=EXP((B+X)/.1417349) 0022 FX=PW-(C I-C2*TWB)/( C3-2. 3*TWB) 0023 FX1=PW/.1417349*((.2646934E-5*TWB-. 12150981E-4*TWJB+.50081415E-2)/X 1 -.1626943E-2) +(C3*C2-2.3*C )/(C3-2.3*TWB)**2 0024 TWB=X1-FX/FX1 0025 IF (A8S(TiB-X1)-.05) 600,4,4 0026 300 ENTH=.24*TDB+HUS*(1035.+. 44*TD8) 0027 GO TO 210 0028 400 TDB=ENTH/(.0058476*ENTH+.290812) 0029 HUS=.00140973*EXP(. 0340657*TDB) 0030 TWB:TDB 0031 500 PVA=BAR0*HUS/f.622+HUS) 0032 600 DENS=(BARO-PVA)/(. 7541*(460.+T0B) ) 0033 1000 RETURN 0034 END TOTAL MEMORY REQUIREMENTS 0005FE BYTES 98

I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 136. 000 + 0- -..+ -.. + _.______ __ + I I I I I * I I I I I I I I I I I I I I I I I I I I I I I * I I I I I I I I I I I I I I I I I *I I I I I I I I I I I I I I I I I I * I I 1* 1 I I I I I * I I I I I I I I I I I I I I I I I *I I I I I I I I I I I I I I I I I I I I I I I I 80. 000 *- - -+-.. -+-_ _._._..+ _._._._.+............+ _..__ _ _ -_ + 0.000. 200.400.600. 800 1. 000 Time (hr) Sample Computer Plot 99

$LI T DRtYDATA i &DATA TOB= i2., I83-=6C., 2 CLdA-=9., 10 IST=85., 3 CTEMP=72., CF-P =22 I., 4 fSP=29-., LVTP=45., 5 H SPMAX=40(0., b I PLO =2, 7 H fMAX =4 0 T M X =60, 8 CONALUC==3. tLVADD=.38 &ENC 9 &DATA CTE,F=72., 1 0 HSP=290., 11 E TP=45., 12 CLOAu=10. &END 1]. &DATA CTEMP=72., 14 HSP=29 0., EVTP=45., i5 CLOAD=6. &END 16 &CATA CTEMF=72., 17 HSP=290., EVTP=45., 1 8 CLO AC= 2. E JNT) 2x L&DATA CTEMP=72., 21 HSP=ZS.,:VTP='45. 22 CLLUAD'=9., C<JIST=7n. = END 23 Ol&UAIA CTE vF=72. 24+ HSP=2cO., L-VTP-45., 2 5 C MOJIST=50*. fLENC 26 &DATA CTEciF=72.,t 27 HSP=2 $., V TP =45., 28 CM!Vul1ST=85., 29 CF M = 25J. AENL) 30 &DATA CTEMP-72., 31 HSP=290., EVTP=45., 32 CFtv= 300. & END 4 &; ATT A T )83=7-2., TBt3b, 41 CLOUAC=9., CMC I ST=5., 4 2 C TE TMP-=72., CF V.=2 20., 43 HSP=290., EV'[P=45., 44 H-SPMAX=40., 45 I PLL T=O, 40 6H AX=45., TBMAX=6/." 4 -+CNADO=3.),K EVAuu=.3L 8 ENC 4d L)DATA CTEMP=72., ~^9 HSP=29g., LEVTP=45. 50 H IAX=35. &tbEND 51 bLATA CTEiF=72., 52 HSP =25'., EVTP=45., 53 H v" X = 4i. 8EN.D ENi J F FILE Sample Data File 100

UNIVERSITY OF MICHIGAN IlEIIH H3D2Um I0I6 3 9015 03525 0698