15radiat.html

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 4.2: Heat Diffusion with Radiation Cooling

Maple Packages for Section 4.2

> restart;

> with( PDEtools ):

>

In Newton's Law of Cooling, it is supposed that if a body with temperature is immersed in a medium with temperature then the body cools according to the model,

T ' = - k ( T - ), T(0) = .

Here, k is the heat transfer coefficient. This model takes the object to be a point source. If we assume that we have a rod that is not insulated, so that we have diffusion and radiation cooling, the equation becomes slightly more complicated. We have an equation with a diffusion term and a radiation cooling term:

We suppose that the ends of the rod are at some constant temperature. For this example, take the end points to have temperature . Also, the rod should have an initial temperature. So as not to be too specific, take the initial temperature to be distributed over the rod as expressed by a function f ( x ). We write the boundary condition and initial conditions as

u( t, 0) = and u( t, L) = , with u( 0, x) = f ( x ).

Here, the length of the rod is L.

This is the problem we address in this section. If were zero, this would be a homogeneous partial differential equation. If is not zero, it is a non-homogeneous partial differential equation .

A review of Section 3.2, Problem 3 might be appropriate here. As indicated in the discussion of that problem, we divide the equation into two parts: the steady state equation and the transient equation. We solve both these equations and find that their sum provides a solution to the original problem. Here are the two equations:

Steady State Equation:

with v(t, 0) = and v(t, L) = .

Transient Equation:

with w(t, 0) = 0 and w(t, L) = 0.

The sum of these two solutions provides a solution for the original problem.

Steady State + Transient

We show here that the steady state solution plus the transient solution provides a solution for the original problem.

> u:=(t,x)->v(x)+w(t,x);

> diff(u(t,x),t);

> diff(v(x),x,x)+k*(v(x)-beta)+subs({diff(w(t,x),t)=diff(w(t,x),x,x)-k*w(t,x)},diff(u(t,x),t));

> collect(%,k);

We see that, because u(t, x) = v(x) + w(t, x), this last becomes

Boundary conditions check, too.

>

We solve these two equations, one at a time.

Solution for the Steady State Equation

The notion of a steady state solution is, first, physical. With a system such as that which modeled radiation cooling, we had a non-homogeneous partial differential equation. The long range forecast was that no matter what the initial condition, the solution for the equation will evolve in time, settling down to a solution that is not changing in time. This solution is steady in time.

Such a notion is easy to characterize mathematically. It must be that the steady solution is, in fact, a solution of the partial differential equation, but also it is not changing in time. The derivative with respect to t is zero.

We say this again with an equation.

0 = ,

with

and .

Since this is a steady state, and time is not a factor, we could consider this as actually an ordinary differential equations in x. We write it with that notation and call solution which is a function of x by v ( x ).

0 = v''(x) - v(x) + k , with s(0) = and s(1) = .

The steady state solution is computed.

> SState:=dsolve({diff(v(x),x,x)-k*(v(x)-beta)=0,v(0)=beta,v(L)=beta},v(x));

We should not be surprised that if the endpoints are held at temperature and the surrounding medium has temperature then the long range forecast is that the rod has temperature .

>

Solution for the Transient Equation

We solve the transient equation. Recall the equation from above:

with w(t, 0) = 0 and w(t, L) = 0.

As in Section 4.1, the method is separation of variables. Suppose that w can be written as

w(t,x) = T( t )X( x ).

In this case, the above equation becomes

T ' X = T X '' - k T X

T '/ T = (X ''- k X) / X.

With the same arguments as used in Section 4.1, this last equation leads to two ordinary differential equations with boundary conditions:

X '' - k X = X, with X(0) = X(L) = 0,

and T ' = T.

We solve these two differential equations. We can change the first one to

X '' = ( k + ) X, with X(0) = X(L) = 0.

From arguments similar to those in Section 3.3, we get the eigenvalues and eigenfunctions for this problem. We find that there are an infinite of possible values:

( k + ) = , n = 1, 2, 3, ...

Corresponding to each n , we have

X ( x ) = sin( n x /L ).

The equation for T becomes

T ' = ( ) T.

Solutions for this equation will be

T ( t ) = exp( - k t ) exp( ).

We check these solutions for the ODE's

> X:=x->sin(n*Pi*x/L);
T:=t->exp(-k*t)*exp(-n^2*Pi^2*t/L^2);

> simplify(diff(X(x),x,x)-k*X(x)+(k+n^2*Pi^2/L^2)*X(x));

> X(0); X(L);

> simplify(diff(T(t),t)+(k+n^2*Pi^2/L^2)*T(t));

>

Having these solutions for the ordinary differential equations, we can construct the general solution for the homogeneous partial differential equation.

u( t, x ) = .

This provides a solution for the transient equation.

We are now ready to make the solution for the original problem. We add the solution for the steady state equation and the equation for the transient equation:

u( t, x) =

In doing our calculations, we use only a finite approximation for this u.

> u:=(t,x)->beta+exp(-k*t)*sum(A[n]*exp(-n^2*Pi^2/L^2*t)*sin(n*Pi*x/L),
n=1..10);

>

Check the solution for the original PDE

We check that this is a solution for the partial differential equation with boundary conditions.

> diff(u(t,x),t)-diff(u(t,x),x,x)+k*(u(t,x)-beta);

> u(t,0); u(t,L);

>

To get the particular solution that satisfies the initial condition

f(x) = u(0, x) = + .

We see that getting the coefficients is simply a matter of computing the Fourier coefficients for the function f(x) - .

At this point, we have not specified a particular L, or k, or , or f. We do this now so that calculations may be made, solutions obtained, and graphs drawn. Take L = 2, = 3, k = 0.005, and

f ( x ) = .

> L:=2; beta:=3; k:=5/1000; f:=x->x^3*(2-x)+3;

> f(x)=u(0,x);

>    for n from 1 to 10 do
   A[n]:=int((f(x)-beta)*sin(n*Pi*x/L),x=0..L)/
            int(sin(n*Pi*x/L)^2,x=0..L);
od;
n:='n':

>

Since we have checked that this u satisfies the partial differential equation and the boundary conditions, we only check that it is an approximation for f ( x ). In fact, we off set the graph of f a little so that we can distinguish these two.

> plot([u(0,x),f(x)+0.01],x=0..L,u=0..5);

>

Finally, we draw a graph of u . What we should expect is that the solution moves from a shape similar to the graph of f to the steady state of .

> plot3d(u(t,x),x=0..2,t=0..1,axes=normal,orientation=[-120,60]);

>

What was new in this section was how to handle a non-homogeneous partial differential equation. The technique was to formulate the steady state equation and the transient equation. It is a good idea to verify that this has been done correctly by checking that the sum of these two is a solution for the original problem. The next step is to solve these two equations, add the solutions, and determine the appropriate coefficients using the Fourier Idea.

Unassisted Maple

We check to see what Maple 8 can do with this equation unassisted. Be aware that we have called in the PDE(tools) at the start of this Section.

> u:='u':

> PDE:=diff(u(t,x),t)-diff(u(t,x),x,x)+k*(u(t,x)-beta);

> ans:=pdsolve(PDE);

$ans := `&where`(u(t,x) = _F1(t)*_F2(x)+3/_C1*(_C1+_C2*exp(1/20*2^(1/2)*x)+_C3*exp(-1/20*2^(1/2)*x)),[{diff(_F1(t),t) = _c[1]*_F1(t), diff(_F2(x),`$`(x,2)) = _c[1]*_F2(x)+1/200*_F2(x)}])$
$ans := `&where`(u(t,x) = _F1(t)*_F2(x)+3/_C1*(_C1+_C2*exp(1/20*2^(1/2)*x)+_C3*exp(-1/20*2^(1/2)*x)),[{diff(_F1(t),t) = _c[1]*_F1(t), diff(_F2(x),`$`(x,2)) = _c[1]*_F2(x)+1/200*_F2(x)}])$

> build(ans);

>

We see that Maple successfully recognizes that the solution should be a sum of a solution for the steady state equation plus a product of functions which will be solutions for the transient equation. Because the solution was obtained without any information about the boundary conditions, this will not be the same solution that we humans got with the assistance of Maple.

Maple handles this problem deftly with other techniques. See Section 8.2.

EMAIL: herod@math.gatech.edu or jherod@tds.net

URL: http://www.math.gatech.edu/~herod