Download Comparison of Five Techniques for Solving Thermal Energy Balance Equations in TRNSYS and more Exercises Differential Equations in PDF only on Docsity!
________________________________________________________________________
CHAPTER
THREE
________________________________________________________________________
SOLVING THE ENERGY BALANCE EQUATION { TC
"SOLVING THE ENERGY BALANCE EQUATION" \l 1 }
After the user specifies which features are to be included in the tank model, the program assembles the energy balance equation. Section 3.1 describes how all the components are assembled into a set of N equations. Five commonly-used techniques for solving simultaneous differential equations are then presented. The chapter provides (or tells where to obtain) the basic equations required to use each technique. Sections 3.2 through 3.5 each describe one particular solution method: providing the key equations and then briefly describing the strengths and weaknesses of the technique. Section 3.6 presents a side-by-side comparison of the five techniques so the user can compare the variation among solution methods.
3.1 Energy Balance Equation{ TC "3.1 Energy Balance Equation" \l 1 } The energy balance equation must account for all possible energy flows into and out of a node. A most general scenario would occur if every energy flow possibility occurred in one single node. Figure 3.1.1 shows all energy flows that could occur in a node.
(U+∆U)As (Tenv-Ti)
Qaux
Node i
UAhx lmtd
² (Ti-1 - Ti)
( k+² k ) X
Ac
UAflue (Tflue - Ti)
² (Ti+1 - Ti)
( k+² k ) X
Ac
m2in Cp (T2in)
or
or
m1in Cp (T1in)
mdown Cp (Ti-1)
mup Cp (Ti)
m1out Cp (Ti)
m2out Cp (Ti)
mup Cp (Ti+1) mdown Cp (Ti) Figure 3.1.1 Energy flows in a node { TC "Figure 3.1.1 Energy flows in a node" \l 1 }
The energy flows designated by the solid arrows represent energy transfer associated with mass flow, and the striped arrows account for energy transfer by conduction, convection or direct energy input. Taking all the energy flows shown in the figure into account, the differential form of the node energy balance is given by equation 3.1.1:
turned on. Hot water to the load is drawn off the top while the cold replacement fluid enters at the bottom. The tank has uniform losses. Using node 2 as an example, equation 3.1.1 reduces to:
M 2 Cp^ dT dt^2 = ( k^^ + ∆^ x∆ 3 → k ) 2 Ac ( T 3 − T 2 )
+ (^ k^ + ∆^ ∆x 1 k →) 2 Ac ( T 1 − T 2 )
+ ( U + 0 ) As ( T env − T 2 )
− γ 2 mÝup Cp ( T 2 )+ γ 4 mÝup Cp ( T 3 )+ γ 5 Qaux 1
Assuming there is flow to the load as shown, γ 2 and γ 4 both have a value of 1, otherwise, they would be zero. If the auxiliary heater was on, γ 5 would be 1, otherwise it would be zero.
Assembling the Equations Once it is known which terms must be included in the energy balance, equation 3.1.1 is computed for each node. Each node's temperature is affected by adjacent nodes, the flue, heaters, heat exchangers, and the environment. With the exception of the node temperatures (and in some cases electric heater power and energy input from the heat exchanger, discussed in later chapters), all the variables in equation 3.1.1 remain constant for the entire TRNSYS time step. Given this, all of the constants can then be grouped together to become a coefficient on each node temperature. Equation 3.1.1 would then take the form of equation 3.1.3:
d Ti d t =^ ai^ Ti −^1 +^ bi^ Ti^ +^ ci^ Ti +^1 +^ di^ (3.1.3)
Rather than storing a, b, c and d as four separate vectors, it is more convenient to store them in a two-dimensional array. The program stores a, b, and c in the matrix A and d in the array d. The computer code (see Appendix A) assembles A and d in the subroutine coeff. For example, the
matrix A would be a banded matrix like the one shown in equation 3.1.4 for a 5 node tank.
dT 1 dt dT 2 dt dT 3 dt dT 4 dt dT 5 dt
b 1 c 1 0 0 0 a 2 b 2 c 2 0 0 0 a 3 b 3 c 3 0 0 0 a 4 b 4 c 4 0 0 0 a 5 b 5
T 1
T 2
T 3
T 4
T 5
d 1 d 2 d 3 d 4 d 5
The program must solve equation 3.1.3 for every node in the tank. The next four sections of this chapter describe solution methods for solving equation 3.1.3, including the strengths and weaknesses of each technique.
3.2 Analytical Solution{ TC "3.2 Analytical Solution" \l 1 } Ideally, equation 3.1.3 would be solved analytically, thus eliminating the inaccuracies associated with using numerical techniques. An analytical solution requires simultaneously solving a set of non-homogeneous differential equations. From equation 3.1.4, let T represent a column vector containing the node temperatures, let T' be the left-hand side of equation 3.1.4, and let d represent a column vector containing d 1 , d 2 , etc. Equation 3.1.4 can then be written in the form of equation 3.2.1:
T' =^ AT +^ d (3.2.1)
If d is ignored for now, equation 3.2.1 is a homogeneous differential equation, the solution to which is:
T ( t^ )^ = e λλ^ t^ ΛΛC (3.2.2)
Where λλ is the vector containing the eigenvalues, (^) ΛΛ is a matrix of eigenvectors, and C is a vector of constants (not to be confused with the lower-case ci , used in equations 3.1.3 and
As will be explained in later chapters, auxiliary heaters may cycle on and off, and, if internal heat exchangers are present, the UA value of the heat exchanger changes with node temperature. Each time something in the energy balance equation changes, a new set of eigenvalues, eigenvectors, particular solutions, and unknown constants must be re-computed. To maintain the goal of a quickly running model, another solution technique is needed. The second reason for not using the analytical method is the mixing algorithm described in section 2.2.4. The mixing of nodes changes their temperatures in a manner that is not included in the energy balance equation. Obtaining an accurate value of (^) Tnew is pointless if the mixing algorithm is going to change the final value. The solution method must work with both T and Tnew , where^ Tnew is the temperature of the node after the solution method^ and^ the mixing algorithm have been performed.
3.3 "DIFFEQ" Solutions{ TC "3.3 "DIFFEQ" Solutions" \l 1 } The old TRNSYS tank model uses the DIFFEQ solution technique to solve the energy balance equations. The solution technique is named as such because it calls the subroutine DIFFEQ, which is part of the TRNSYS software. DIFFEQ returns an analytical solution to equation 3.3.1.
dTi dt =^ ai^ Ti^ +^ bi^ (3.3.1)
Where (^) ai would be the same as (^) bi in equation 3.1.3 and (^) bi would be
bi = a 2 i ( T (^) i − 1 + Tnew,i − 1 )+ c 2 i ( T (^) i + 1 + Tnew,i + 1 )+ di (3.3.2)
The coefficient (^) bi is not constant, but is assumed so in the solution. The analytical solution to equation 3.3.1 for a constant value of bi is:
Tnew,i = Ti (t) (^) + b aii eai∆t^ − b aii (3.3.3)
The average temperature over the time step t is defined as:
Tavg,i = (^) ∆^1 t^ Ti (t ) + b aii eait^ − b aii dt 0
∆t ∫ (3.3.4)
Integration of equation 3.3.4 yields:
Tavg,i =
Ti ( ) t +^ b aii
ai ∆t ( e^ bi^ ∆t^ −^1 )−^
bi ai (3.3.5)
For every time step t, the new node temperatures Tnew,i are calculated by assuming all adjacent nodes (nodes (^) i − (^) 1 and (^) i + (^) 1) are at a time-averaged temperature. The program obtains the average temperatures by iterating through the nodes (holding time constant) until the change in Tnew,i is within a specified tolerance. Although iterating with average node temperatures works, using time-averaged adjacent node temperatures in equations 3.3.3. and 3.3.5 to solve the energy balance may be an incorrect assumption. The problem lies in using temperatures in the same equation that occur at different times. Time-averaged adjacent node temperatures are computed and then become part of the constant bi in equation 3.3.2. However, equations 3.3.3 through 3.3.5 are based on the initial temperature of the node i. Consider, for example, flow coming into a tank from below. The flow into node i would be at the time-averaged temperature of the node below it. However, the flow leaving node i (going into the node above) would be at its initial temperature, because equations 3.3.3 through 3.3.5 do not permit the use of the time-averaged temperature for the node being solved for. Although using the average temperature for adjacent nodes is not entirely correct, this technique does seem to work. Still, one may wish to proceed with caution when using the DIFFEQ method. The advantage of this method is its ability to take large time steps and still
when compared against other solution techniques on an equal time step size basis, forward Euler is always less accurate.
3.5 Crank-Nicolson Solution{ TC "3.5 Crank-Nicolson Solution" \l 1 } For reasons that will become evident in the next section, the Crank-Nicolson solution was found to be the best overall technique for solving the energy balance equations. Crank- Nicolson, like DIFFEQ, requires iterations to obtain an average temperature, but can consequently use large time steps. Unlike the Euler solution, it is unconditionally stable. The new node temperatures are computed using equation 3.5.1.
Tnew,i = ∆t
ai 2 ( T^ i −^1 +^ Tnew,i^ −^1 )+^
bi 2 ( T^ i^ +^ Tnew,i )+ ci 2 ( T^ i^ +^1 +^ Tnew,i +^1 )+^ di
+ Ti (3.5.1)
The computer loops through equation 3.5.1 until the change in Tnew,i is less than some specified value. Depending on the particular situation, 2-3 iterations are usually required to converge to within 0.001 C. The advantage to the Crank-Nicolson solution is its ability to take large time steps and still obtain an accurate solution. The problem lies with the average temperature, which is computed using:
Tavg = ( T^ new,i 2^ +^ Ti ) (3.5.2)
As mentioned earlier, (section 2.2.4) the average temperature is the key element for accurate tank energy balance calculations. The problems incurred by using inaccurate average node temperatures will be addressed in chapter 6. When the Crank-Nicolson technique is compared against the Euler technique, the time step can be much larger. However, regardless of solution
method used, increasing the size of the time step increases the error in both the average and final temperatures.
3.6 Solution Method Comparison{ TC "3.6 Solution Method Comparison" \l 1 } Using a three-node tank as a basis, the analytical, DIFFEQ, DIFFEQ2, forward Euler and Crank-Nicolson solution methods were compared. Assume the energy balance equation for this tank is:
dT 1 dt = −^1 T^1 +^2 T^2 +^0 T^3 +^1 dT 2 dt =^2 T^1 −^1 T^2 +^1 T^3 +^2 dT 3 dt =^0 T^1 +^1 T^2 −^0_._^5 T^3 +^10
The critical Euler time step is:
∆tcrit = minimum − b^1 i
i = 1 ,n
= (− − 11 ) = 1 hour (3.6.2)
Where bi is the constant coefficient on Ti. (see section 3.1). For accuracy reasons, it may not be desirable to use a time step that is greater than 1/6 the critical Euler value (see chapter 4). For the comparison described below, the value after 1 hour was examined. The tank had a temperature of 0 at time=0. The coefficients in equation 3.6.1 were chosen to cause the bottom node to be warmer than the other nodes. This causes the need to include the mixing algorithm (see section 2.2.4). The mixing of the nodes causes the value of Tnew,i to differ from the equation's solution; and mixing of nodes also makes it impossible to solve the energy balance in the true analytical fashion described in section 3.2, because mixing de-couples the methods by which energy may enter (or leave) a node. Tnew,i is the temperature at the end of the time step
iteration. The iterations / time step column shows on average how many times the computer must loop through the equations to get the change in (^) Tnew less than 0.0001 C. The solution method chosen for use in the new tank model was the Crank-Nicolson. The computer code used to calculate the numbers in table 3.6.1 is listed in Appendix C.
Comparison Without Mixing Algorithm The mixing algorithm changes the final result of solution techniques that might otherwise work fine. Equation 3.6.1 has been solved using the five methods without the mixing algorithm included. Table 3.6.2 presents the results. In this case, not all nodes will be at the same temperature. For the sake of brevity, only the results for the top and bottom nodes are presented in the table.
Table 3.6.2a Solution method comparison without mixing algorithm, top node { TC "Table 3.6.2a Solution method comparison without mixing algorithm, top node" \l 1 }
time step=1/6 Solution Iterations / %error obtained time step METHOD Theoretical 4.8177 0. DIFFEQ 4.8459 7.0 0. DIFFEQ 2 3.0629 1.0 36. Euler 3.5616 1.0 26. Crank-Nicolson 4.8678 6.5 1.
time step=1/60 Solution Iterations / % error obtained time step METHOD Theoretical 4.8177 0. DIFFEQ 4.8175 3.0 0. DIFFEQ 2 4.5926 1.0 4. Euler 4.6706 1.0 3. Crank-Nicolson 4.8177 3.0 0.
As can be seen in the tables, the DIFFEQ method is slightly more accurate at larger time steps, but at the expense of a few more required iterations. When the size of the time step decreases, the Crank-Nicolson solution method becomes more accurate than the DIFFEQ method.
Table 3.6.2b Solution method comparison without mixing algorithm, bottom node { TC "Table 3.6.2b Solution method comparison without mixing algorithm, bottom node" \l 1 }
time step=1/6 Solution Iterations / %error obtained time step METHOD Theoretical 10.2348 0. DIFFEQ 10.2584 7.0 0. DIFFEQ 2 9.2335 1.0 9. Euler 9.7089 1.0 5. Crank-Nicolson 10.2708 6.5 0.
time step=1/60 Solution Iterations / % error obtained time step METHOD Theoretical 10.2348 0. DIFFEQ 10.2345 3.0 0. DIFFEQ 2 10.1058 1.0 1. Euler 10.1669 1.0 0. Crank-Nicolson 10.2346 3.0 0.
