Objectives:
- Math: Time dependent heat transfer in 2D
- Computing: Animation
- Science: Solar energy storage
General Information: This lecture uses a modification of the previous SOR procedure to approximate the solution of the time dependent heat equation in 2D. We use the implicit time discretization.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-4-94
Copyright ©1994, R. E. White. All rights reserved.
3.6 Time Dependent Problems
3.6.1. Introduction.
In this lecture we consider a time dependent 2D heat diffusion problem. We will use the implicit time discretization method which will generate a sequence of problems that are similar to the steady state 2D heat equation. Here we will use a modified version of the SOR procedure to approximate the solutions of this sequence.
Consider a passive solar storage unit which collects energy by the day and gives it back at night. A simple example is a thick concrete floor in front of windows which face the sun during the day. Suppose it is relatively long, has depth H and width L. So, we can neglect diffusion in the long direction and assume diffusion in the y direction (H) and the x direction (L). Since the surrounding temperature will vary with time, the amount of heat that diffuses in and out of the top will depend on the time. The problem is to determine the effectiveness of the passive unit as a function of its geometric and thermal properties.
 |
| Figure: Cross Section of Passive Solar Storage Unit |
The model has the form of a time and space dependent partial differential equation. The diffusion is again governed by the empirical Fourier heat law. For the non steady state problem we must consider the change in the amount of heat energy that a mass can retain as a function of temperature. If the temperature varies over a large range, or if there is a change in physical phase, then this relationship is nonlinear. However, for small variations in the temperature we have
- change in heat =
cp(u - u ) volume where
is the density,
cp is the specific heat and
u - u is the change in temperature.
When this is coupled with diffusion in both the x and y directions for the volume (dx dy T), we get
- cp(u - u ) (dx dy T) =
f(x, y, t+dt) (dx dy T) dt
+ (dx T) dt (Kuy(x, y+dy, t+dt) - Kuy(x, y, t+dt))
+ (dy T) dt (Kux(x+dx, y, t+dt) - Kux(x, y, t+dt)).
This approximation gets more accurate as dx, dy and dt go to zero. So, divide
by
Time Dependent 2D Heat Diffusion.
- cp
ut = f(x,y) + (Kux(x,y))x +
(Kuy(x,y))y
(1) - for (x,y) in (0,L)×(0,H),
- f(x,y) is the internal heat source,
K is the thermal conductivity.Boundary conditions: Kux = 0 for x = L, u = 60 for x= 0,
-Kuy = 0 for y = 0 and
Kuy = c(usur(t) - u) for y = H with
usur(t) = 60 + 30 sin(t Pi/12).(2) Initial condition: u = 60 for t = 0. (3) - for (x,y) in (0,L)×(0,H),
The derivation of (1) suggests the implicit time discretization method. Let
k denote the time step with
-
cp(uk+1 - uk) (dx dy T) =
f(x,y,(k+1)dt) (dx dy T) dt + (dx T) dt (Kuyk+1(x,y + dy) - Kuyk+1(x,y)) + (dy T) dt (Kuxk+1(x+dx, y) - Kuxk+1(x,y)).
The space variables can be discretized just as in the steady state heat diffusion problem. Thus, for each time step we must solve a linear algebraic system where the right side changes from one time step to the next and the right side equals
Implicit Time Discretization Algorithm.
In the calculations below we have used a variation of the SOR algorithm to
solve each linear algebraic system. This is done in the procedure BVP2D where
the boundary conditions with the derivatives have been inserted in the SOR
loop given by the loop index m.
The output data is given by an animation using the Maple procedure display
with the graphic structures give by the Maple procedure matrixplot. The
initial temperature was uniformly set at 60, and four of the 24 frames are in
the figures for times 6, 12, 18 and 24. In these figure the x axis is along
the rows, and the front right in the figure is the top of the storage unit.
The y axis is along the columns, and the back right in the figures has a fixed
temperature equal to 60. Heat is entering through the top at times 6 and 12.
Heat is leaving through the top at times 18 and 24.
Maple Code for Passive Solar Energy Storage.
The effectiveness of the storage unit must be a function of the heat stored,
the cost of construction and the cost of alternative heat sources.
The amount of heat stored will depend on the geometric
and thermal parameters of the storage unit as well as the source of heat from
the sun. All these parameters can be varied in this code, and this will allow
for a systematic investigation of many different configurations without having
to construct the storage unit. Of course, the storage unit is just one
segment of the entire structure, and one should take other parts of the
structure into consideration.
The implicit time discretization is called implicit because one must solve an
algebraic linear system to find the temperature at the next time step. This
contrasts with the explicit time discretization where the temperature at the
next time step can be found without solving an algebraic system. The advantage
of the implicit method is that there is no stability condition on the time
step.
Explicit Time Discretization Algorithm.
cpuk/dt +
f(x,y,(k+1)dt).
uold = u(x,y,0)
for k = 1, maxk
approximate the solution by the finite difference method:
cp (u - uold)/dt = f(x,y) + (Kux(x,y))x + (Kuy(x,y))y
with appropriate boundary conditions
uold = u
endloop.
Input Data:
> with(linalg):
> with(plots):
> L:=10:
> Time:=24:
> nt:=24:
> rho:=119:
> cht:=.21:
> cond:=.81:
> eps:=.1:
> nx:= 10:
> ny:=6:
> c:=2:
> H:=6:
> w:=1.85:
> usur:=(t)->60 +30* sin(Pi*t/12);
Procedures BVP2DTIME and BVP2D are Defined:
> BVP2DTIME:=proc(u,uold,c,usur,Time,nt,rho,cht,cond,L,nx,ny,H,w,eps):
> u:=matrix(nx+1,ny+1,60):
> uold:=matrix(nx+1,ny+1,60):
> uu.0:=matrixplot(uold,style=wireframe,axes=boxed):
> for k from 1 to nt do
> BVP2D(u,uold,c,usur,k,Time,nt,rho,cht,cond,L,nx,ny,H,w,eps):
> uu.k:=matrixplot(u,style=wireframe,axes=boxed):
> uold:=u:
> od:
> end:
> BVP2D:=proc(u,uold,c,usur,k,Time,nt,rho,cht,cond,L,nx,ny,H,w,eps):
> h:=H/ny:
> dt:= Time/nt:
> maxit:=400:
> tol:=eps*h*h:
> cc:=rho*cht/dt:
> for m from 1 to maxit do
> numi:=0:
> j:= 1:
> for i from 2 to nx do
> utemp:= evalf((uold[i,j]*h*h*cc + cond*(2*u[i,j+1]
> + u[i+1,j] + u[i-1,j])))/(4*cond+h*h*cc):
> utemp:= (1-w)*u[i,j] + w*utemp:
> error:= abs(utemp - u[i,j]):
> u[i,j]:=utemp:
> if error < tol then numi:=numi +1
> fi:
> od:
> i:=nx+1:
> utemp:= evalf((uold[i,j]*h*h*cc + 2*cond*u[i,j+1]
> + 2*cond*u[i-1,j] ))/(4*cond+h*h*cc):
> utemp:= (1-w)*u[i,j] + w*utemp:
> error:= abs(utemp - u[i,j]):
> u[i,j]:=utemp:
> if error < tol then numi:=numi+1:
> fi:
> for j from 2 to ny do
> for i from 2 to nx do
> utemp:= evalf((uold[i,j]*h*h*cc + cond*(u[i,j-1] + u[i-1,j]+
> u[i+1,j] +u[i,j+1] ))/(4*cond+h*h*cc)):
> utemp:= (1-w)*u[i,j] + w*utemp:
> error:= abs(utemp - u[i,j]):
> u[i,j]:=utemp:
> if error < tol then numi:=numi + 1
> fi:
> od:
> i:=nx+1:
> utemp:= evalf((uold[i,j]*h*h*cc+ cond*(u[i,j-1]
> +u[i,j+1]+2*u[i-1,j] ))/(4*cond+h*h*cc)):
> utemp:= (1-w)*u[i,j] + w*utemp:
> error:= abs(utemp - u[i,j]):
> u[i,j]:=utemp:
> if error < tol then numi:=numi+1:
> fi:
> od:
> j:= ny + 1:
> for i from 2 to nx do
> utemp:= evalf((uold[i,j]*h*h*cc+2*c*h*usur(k*dt)
> +cond*(2*u[i,j-1] + u[i-1,j] + u[i+1,j]))/
> (4*cond +h*h*cc + 2*c*h)):
> utemp:= (1-w)*u[i,j] + w*utemp:
> error:= abs(utemp - u[i,j]):
> u[i,j]:=utemp:
> if error < tol then
numi:=numi+1:
fi:
> od:
> i:=nx+1:
> utemp:= evalf((uold[i,j]*h*h*cc+2*c*h*usur(k*dt)
> +cond*(2*u[i,j-1] + 2*u[i-1,j]))/
> (4*cond +h*h*cc + 2*c*h)):
> utemp:= (1-w)*u[i,j] + w*utemp:
> error:= abs(utemp - u[i,j]):
> u[i,j]:=utemp:
> if error < tol then numi:=numi+1:
> fi:
> if numi=(nx)*(ny+1) then break fi:
> od:
> end:
BVP2DTIME is Called and the Output Data is Animated:
> BVP2DTIME(u,uold,c,usur,Time,nt,rho,cht,cond,L,nx,ny,H,w,eps):
> LL:=[seq(uu.l,l=0..nt)]:
> display(LL,insequence=true);
Figure: Temperature at Time = 6
Figure: Temperature at Time = 12
Figure: Temperature at Time = 18
Figure: Temperature at Time = 24
uold = u(x,y,0)
for k = 1, maxk
approximate the solution by the finite difference method:
cp (u - uold)/dt = f(x,y) + (Kuoldx(x,y))x + (Kuoldy(x,y))y
with appropriate boundary conditions
uold = u
endloop.