Ch9.3 Model for a Hot Water Heater

• Let \( U(t) \) = temperature of the water at time t.
• Let \( u_0 \) = initial temperature of the water,
• \( u_f \) = final temperature,
• \( m \) = mass of water in the heater
• \( q \) = rate of heat energy supplied
• \( S \) = surface area \( m^2 \) from which heat can escape, which in the case of the hot water heater is the surface of the tank.

• Assume the water in tank is well stirred so that temperature remains homogeneous throughout.
• Otherwise temperature would be function of time and position
• We would need to account for heat transfer (by conduction and convection) between different points in tank.
• By assuming that the water is well stirred, we only have to consider the temperature of the water as a function of time.

• Start with input-output diagram for the heat content of the system by applying the balance law.

\[ \begin{Bmatrix} \mathrm{rate \, of \, change} \\ \mathrm{of \, heat } \end{Bmatrix} = \begin{Bmatrix} \mathrm{rate \, heat \, produced} \\ \mathrm{by\, heating \, element} \end{Bmatrix} - \begin{Bmatrix} \mathrm{rate \, heat \, lost} \\ \mathrm{to \, surroundings} \end{Bmatrix} \]

• Net heat produced goes into raising temperature of system.
• Need mathematical expression for heat required to change temperature.
• Use fundamental equation relating heat to temperature via specific heat (Ch9.2):

\[ Q = cm \frac{dU}{dt} \]

• Here \( Q \) is rate of change of heat content, \( c \) is specific heat of water, \( m \) is mass of water and \( U(t) \) is temperature at time t.

• On the previous slide, we found an expression for left side of word equation:

\[ Q = cm \frac{dU}{dt} \]

• We next determine expressions for the right side.

\[ \begin{Bmatrix} \mathrm{rate \, of \, change} \\ \mathrm{of \, heat } \end{Bmatrix} = \begin{Bmatrix} \mathrm{rate \, heat \, produced} \\ \mathrm{by\, heating \, element} \end{Bmatrix} - \begin{Bmatrix} \mathrm{rate \, heat \, lost} \\ \mathrm{to \, surroundings} \end{Bmatrix} \]

• Assume heating element produces heat at a constant rate per unit time \( q \):

\[ \begin{Bmatrix} \mathrm{rate \, heat \, produced} \\ \mathrm{by\, heating \, element} \end{Bmatrix} = q \]

• Word equation from before:

\[ \begin{Bmatrix} \mathrm{rate \, of \, change} \\ \mathrm{of \, heat } \end{Bmatrix} = \begin{Bmatrix} \mathrm{rate \, heat \, produced} \\ \mathrm{by\, heating \, element} \end{Bmatrix} - \begin{Bmatrix} \mathrm{rate \, heat \, lost} \\ \mathrm{to \, surroundings} \end{Bmatrix} \]

• Rate of heat lost to surroundings: Newton's law of cooling

\[ \begin{align*} \begin{Bmatrix} \mathrm{rate \, heat \, lost} \\ \mathrm{to \, surroundings} \end{Bmatrix} & = hS \begin{Bmatrix} \mathrm{temperature} \\ \mathrm{difference} \end{Bmatrix} \\ \\ & = hS\left[ U(t) - u_s \right] \end{align*} \]

• Here \( S \) is the surface area of heater, \( h \) is the Newton cooling coefficient, \( U(t) \) is the water temperature and \( u_s \) is the temperature of surroundings.

\[ \begin{Bmatrix} \mathrm{change \, in} \\ \mathrm{heat \, energy } \end{Bmatrix} = \begin{Bmatrix} \mathrm{amount \, heat } \\ \mathrm{produced \, by } \\ \mathrm{heating \, element} \end{Bmatrix} - \begin{Bmatrix} \mathrm{amount\, heat} \\ \mathrm{lost \, to } \\ \mathrm{surroundings} \end{Bmatrix} \]

• Change in heat energy = \( mc \left[ U(t + \Delta t) - U(t) \right] \)
• Amount heat produced = \( q \Delta t \)
• Amount heat lost to surroundings = \( hS \left[ U(t^{*}) - u_s \right]\Delta t \)
• \( U(t^{*}) \) = temperature at \( t^{*} \in [t, t + \Delta t] \)

\[ \begin{Bmatrix} \mathrm{change \, in} \\ \mathrm{heat \, energy } \end{Bmatrix} = \begin{Bmatrix} \mathrm{amount \, heat } \\ \mathrm{produced \, by } \\ \mathrm{heating \, element} \end{Bmatrix} - \begin{Bmatrix} \mathrm{amount\, heat} \\ \mathrm{lost \, to } \\ \mathrm{surroundings} \end{Bmatrix} \]

\[ \begin{align*} mc \left[ U(t + \Delta t) - U(t) \right] &= q \Delta t - hS \left[ U(t^*) - u_s \right]\Delta t \\ mc \lim_{\Delta t \rightarrow 0} \frac{U(t + \Delta t) - U(t)}{\Delta t}\Delta t &= q - hS \lim_{\Delta t \rightarrow 0} \left[ U(t^*) - u_s \right] \\ cm \frac{dU}{dt} &= q - hS\left[ U(t) - u_s \right] \end{align*} \]

\[ \frac{dU}{dt} = \frac{q}{cm} - \frac{hS}{cm}\left(U - u_s \right) \]

• Units for left side of ODE are clearly \( \frac{C}{sec} \).
• Units for right side ODE terms (see previous slide also):

\[ \begin{aligned} \left[\frac{q}{cm}\right] & = \frac{ \frac{J}{sec} } { \left(\frac{J}{kg ~C}\right)kg } = \frac{C}{sec} \\ \left[\frac{hS}{cm}\left(U - u_s \right)\right] & = \frac{ \left(~\frac{J}{m^2 sec~ C}\right)\left(m^2~\right)\left(~C~\right) } { \left(\frac{J}{kg ~C~}\right)~kg~ } = \frac{C}{sec} \end{aligned} \]

Ch93Ex1 <- function(T) {
#Ch9.3: Perform Rk4 for water heater
#T is time length in hours for [0, T]
N <- 3600*10   #N is number of time nodes
h <- 3600*T/N    #Time step size in seconds
#System Parameters
t0 <- 0
u0 <- 15 
us <- 15
q <- 3600    
c <- 4200    
m <- 250      
h2 <- 12  #h is used above for step size   
S <- 3.06    
#Slope function from ODE
f <- function(U) {q/(c*m) - (h2*S)/(c*m)*(U - us)}  

#Initialize vectors for time t and temperature U. 
t <- rep(0, N)
U <- rep(0, N)
t[1] <- t0
U[1] <- u0
#Runge-Kutta Loop (Generate temperature vector U)
for(i in 1:N) {
a1 <- h*f(U[i]);          #f = slope of U
b1 <- h*f(U[i]+0.5*a1);   #Half-step predictor 
c1 <- h*f(U[i]+0.5*b1);   #Half-step predictor 
d1 <- h*f(U[i]+c1);       #Full-step predictor 
U[i+1]<-U[i]+(1/6)*(a1+2*b1+2*c1+d1);#Next U value
t[i+1] <- t[i] + h  #Next t value, need for plot}