In this exercise you will solve numerically the 1D Heat Equa

In this exercise, you will solve numerically the 1-D Heat Equation. Consider a beam that is heated up in the center: Beam heated up in the centre Uneven Temperature distribution along the beam At time t 0, the heating stops. At that time, the temperature distribution along the beam is determined as: which will be the initial condition for your computation (with To-100K and Tbase 300K). As time elapses, the heat from the centre of the beam will be conducted towards the ends. The ends of the beams are cooled convectively with a convective heat transfer coefficient of h-85 W/m2 K, and an external temperature of T-300K, which determines the boundary condition for your computation. The heat diffusivity is a 0.5. You will now use the 1-D Heat Equation to monitor the progress of the temperature distribution. There is no convective or radiative heat transfer along the beam Use the Finite-Difference discretization from Exercise 5. Write a program in a programming language of your choice, such as C++, C, Fortran or Python. You can use Matlab, but only if you do not use any of the higher order functions to solve PDEs.

Solution

function advancedifference=heattransfer(T,dydt)

%our principal equation is T(x)=T0*sin(x*pi()/L)+Tbase

%for this problem we need to find all the differential equation:

%for a convective heat transfer, we have that the heat transfered is:

%Q=m*h*dT/dt, rewritting in terms only of temperature we have q=h*dT/dt

%so our differential equation will be h*dT/dt=T0*sin(x*pi/L)/dt+Tbase/dt

%we are using advanced finite differences to solve this, because we have

%not information about the finite differences used in \"Excersice 5\".

%To determine boundary conditions we need to solve the first equation we

%have for x=0 and x=L. The function sin is 0 in x=0 and 0 in x=L.

%So we have the boundary conditions are Tbase at both sides.

%then, to make finite difference we have to stablish a value for L that

%will not affect our results, lets assume L=20.

T0=100;

Tbase=300;

F=inline(T);

x0=0;

h=20/20;

b=20;

a=0;

n=(b-a)/h;

%this give us the rate of cooling of every node

for i=1:n

advancedifference(i)=(F(i+1)-F(i-1)/(2*h));

x0=x0+h;

end

advancedifference

%for L/2 we have the following equation

t=0.5;

F=inline(dydt);

a=0;

b=1;

n=(b-a)/t;

x=a:h:b;

y0=300;

t0=0;

for i=1:n

  

k1=F(t0);

k2=F(t0+0.5*h);

k3=F(t0+0.5*h);

k4=F(t0+h);

y(i+1)=y0+(1/6)*(k1+2*k2+2*k3+k4)*h;

t0=t0+h;

end

y;

plot(y);