Time dependent heat transfer

18 views (last 30 days)
Bruce Griffin on 29 Sep 2021
Answered: William Rose on 30 Sep 2021
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)
%pause(1)
end
Tn=Tn.*400./pi;
figure (2)
surf(X,Time,Tn)
xlabel('x (cm)')
ylabel('time (s)')
zlabel('Temperature')

William Rose on 30 Sep 2021
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 on 30 Sep 2021
This is exaclty what I was looking for. Thank you very much.

William Rose on 30 Sep 2021
Edited: William Rose on 30 Sep 2021
@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.

William Rose on 30 Sep 2021
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 on 30 Sep 2021
@Bruce Griffin, you're welcome.

William Rose on 30 Sep 2021
@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.

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!