Time dependent heat transfer

Im trying to code and simulate the following problem.

The two sides are suppose to be zero and the third is 100 degrees

However when I run my code my temperature is becoming way to high on scales of x10^4. Is there something I am missing?

clear all

xmax=10;

timemax=20;

numx=1000;

numtime=1000;

x=0:xmax/numx:xmax;

time=0:timemax/numtime:timemax;

[X,Time]=meshgrid(x,time);

alpha=.7;

nmax=1;

Tn=zeros();

for n=1:2:nmax

Tnew=(1/n)*exp(((n*pi*alpha/10)^2)*-Time)*sin(n*pi*X/10);

Tn=Tn+Tnew;

%figure(1)

% surf(x,time,Tn)

% view(0,90)

%shading interp

%pause(1)

end

Tn=Tn.*400./pi;

figure (2)

surf(X,Time,Tn)

shading interp

xlabel('x (cm)')

ylabel('time (s)')

zlabel('Temperature')

Answer 1

[William Rose]

@Bruce Griffin,

Change nmax from 1 to 500.

I really like how you used meshgrid to generate all the time, position combintions as arrays Time and X, and this allows you to compute the array Tnew without using a for loop for time or a for loop for position. Nicely done!

Change

Tnew=(1/n)*exp(((n*pi*alpha/10)^2)*-Time)*sin(n*pi*X/10);

to

Tnew=(1/n)*exp(((n*pi*alpha/10)^2)*-Time).*sin(n*pi*X/10);

The key was to replace matrix multiplication, "*", with element-wise multiplication, ".*" .

The plot is good now:

It is 0 degrees at x=0 and x=10, as desired. It is 100 across the bar at time t=0. It loses heat gradually with time. The temperature distribution across the bar becomes mroe rounded with time. All these results are what we expect. If you look closely at the temperature at time 0, you see "ringing" at the edges x=0 and x=10. This is the expected Gibbs phenomenon, which arises when a square wave is represented by a finite Fourier sum. If you make nmax smaller, the Gibbs phenomenon becomes more obvious.

All of this means that now you have code that computes the answer equation. But this is not a simulation. For a simulation, you should write another script, or add to this script. The new code should solve the difference equations which I derived in my previous post, for each successive time step, given the initial condition and given the boundary conditions. In this way you will fill in a surface like the surface above, starting at t=0 and progressing by with each step.

In my previous answer, I assumed that i was the index for position and j was the index for time. But your Tn is the other way around. So adjust accordingly.

[Bruce Griffin]

This is exaclty what I was looking for. Thank you very much.

Answer 2

[William Rose]

@Bruce Griffin, You said the third [side] is at 100. That is not true. This problem has one spatial dimension and one time dimension. The 10 cm long bar is at a uniform temperature of 100 initially. The ends are clamped to temperature 0 at time 0.

You said you are supposed to simulate the system. The code you have posted is not a simulation. It is just an attempt to compute the solution using the answer equation. I'm not sure exactly why your code gives bad answers. One obvious problem is that the loop for i=1:2:nmax executes one time only, because nmax=1.

A simulation is a numerical solution to the partial differential equation for heat. What is the partial differential equation for heat?

where u(x,t) is the temperature at position x and time t. How does it translate to a discrete system?

therefore

Now let i be the index for position, let j be the index for time. Then

Use the equation above to simulate the time evolution of the temperature distribution.

See if it matches the equation you have been given.

Answer 3

[William Rose]

@Bruce Griffin,

I added code to do an actual simulation. I found an interesting and instructive result. When numx=1000 and numtime=1000 (your values), the first order partial difference equaiton is unstable. In other words, simulated temperature Ts(i,j) starts to oscillate, and the oscillations grow bigger with time. This was evident by the values of -NaN and +NaN in array Ts, which caused errors when I tried to plot Ts. The solution is to use a larger , a smaller , or both. So I changed to numx=500, numtime=2000. Still unstable. I continued until I got to stability: numx=100, numtime=10000. This resulted in a surface plot of Ts that was qualitatively similar to the plot of Tn, the answer equation. See first plot. However, the quantitative agreement between Ts and Tn is not great. See the second plot, with red and blue lines.

The temperature across the bar is declining more quickly in the simulation than in the answer equation. Maybe a larger numx would help, but if I raise numx, I will probably need to raise numtime to maintain numerical stability. This will cause array Ts(i,j) to be quite large, and could cause out-of-memory errors or very long run times.

[William Rose]

The answer provided in the book (Screen Shot 2021-09-29 at 2.37.05 PM.png in the original post) is wrong, and the simulation is right! That is why the simulation results are different from the answer equation. When the answer equation is corrected, the answer equation surface changes, and matches the simulation surface.

Dimensional analysis shows that the top equation is wrong. "10" in the equations above has units of length, t has units of time, thermal diffusivity α has units of length squared over time, and n and π are dimensionless. The power that e is raised to, in the top equation, has units of , and it should be dimensionless. The power that e is raised to in the lower equation is dimensionless, as required. Also, see https://tinyurl.com/4kfbazyb. The solution to this problem is given in dimensionless units, in equation 25 on page 8. The equations for de-dimensionalizaton on the earlier pages, if applied in reverse to equation 25, yield the correct equation above.

[Bruce Griffin]

I really appreciate your answer and help. Yeah I guess I missed represented what I was trying to accomplish. However I took your suggestions and expanded my answer. So thank you

[William Rose]

@Bruce Griffin, you're welcome.

Answer 4

[William Rose]

@Bruce Griffin, here is the code that does the simulation. The original incorrect answer equation is in the code, but it is commented out. YOu can use commenting to select the correct or incorrect equation, and you will see that the "answer" and the simulation match for the correct equation and not for the other.

The fact that simulation reveals a subtle error in the "answer" is a very nice demonstration of the importance of simulation. The error in this case is subtle because α is close to unity (). Therefore there's not a huge difference between having α or in the exponent. If α had been set to unity, we never would have realized there was an error in the orignal answer equation.