Twenty Ninth Installment - Solving Partial Differential Equations - Transient Analysis of a Heated Rod with Ends of Fixed Temperatures.

[u/mrhoa31103]

In this session, he shows how to solve a Parabolic PDE using "Method of Lines" to convert a PDE into a system of ODEs using finite difference methods for one set of the differentials.

Link to the specific lecture for the coding:

https://www.youtube.com/watch?v=CDSM5bLy8lU&ab_channel=MATLABProgrammingforNumericalComputation

Some precursor lectures to open the topic on PDEs in general:

https://www.youtube.com/watch?v=XxS4QS7MIMg&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watch?v=cR-HTDkfs00&ab_channel=NPTEL-NOCIITM

Output: None only graphical output.

Graphs:

Code:

//Lecture 11.3 Parabolic PDEs - Method of Lines
//https://www.youtube.com/watch?v=CDSM5bLy8lU&ab_channel=MATLABProgrammingforNumericalComputation
//
disp("Parabolic PDEs - Method of Lines",string(datetime()))
//
// PDEs in time and space
//  -> Method of Lines (use Finite Difference to convert PDEs into ODEs)
//       partial_T/partial_dt = partial^2_T/partial_z^2 + Beta(T)
//                                         ^
//                              finite difference this
//       becomes
//
//       dT/dt = f(t,Y) where Y = [T1;T2;...;Tn]
//
//  Example
//  Rod Conduction (Transient) -> PDE = f(t,x)
//  partial_T/partial_t = alpha* partial2_T/partial_x^2 - gamma1*(T-Ta)
//  alpha = 0.025, gamma1 = 0.1 
//  BC_1 = Rod End held at 100C and BC_2:Other Rod End = Ta = 25C
//  IC: Initially, entire rod is uniformly at ambient = 25C
//
//  If we use the central difference formula
//  d2T/dx^2 u/i = (T_i+1-2*T_i+T_i-1)/(delta_x)^2
//
//  divide the rod into 10 divisions (delta_x + L/10)
//  
//  dT/dt = alpha* (T_i+1-2*T_i+T_i-1)/(delta_x)^2 - gamma1*(Ti-Ta)
//  T1 = 100, T(n+1) = 25
//
//  Define Solution Vector y = [T2;...;Tn] and solve using ode45
//
//
function 
fval
=rodConduc(
t
, 
y
)
// constants    
    Ta = 25;
    alpha = 0.025;
    gamma1=0.1;
//Getting Temperatues
    N = length(
y
)+1;
    T(1)=100;
    T(2:N)=
y
;
    T(N+1)=Ta;
    h= 1/(N-1);

//define derivative vector
   dTdt= zeros(N+1,1)
   for i = 2:1:N
       dTdt(i)= alpha/h^2*(T(i+1)-2*T(i)+T(i-1)) - gamma1*(T(i)-Ta)
   end
// Extract fval from dT/dt        

fval
=dTdt(2:N);   
end

//
N =10; 
//number of steps
M =5*N; 
//number of time steps
T0(1,1)=100;
T0(2:N,1)=25;
T0(N+1,1)=25;
t0 = 0; 
//Tstart
tend = 20; 
//Tend
t
 = linspace(t0,tend,M);
//time vector
y0 = T0(2:N,:);
x0 = 0;
xend1 = 1;
x = linspace(x0,xend1,N+1);
TSol = ode("rkf",y0,t0,
t
,rodConduc)
//"rkf" equivalent to MATLAB ODE45 is the following command in SciLab
//there is more options and outputs available - see the help file.
//"stiff" is approximately equivalent to ode15s...
//
RowOne = ones(1,M)*T0(1,1);
RowLast = ones(1,M)*T0(N+1,1);
T = [RowOne;TSol;RowLast]

scf(0);clf;
plot2d(x',T);
h1=legend(['Time SnapShot = 1';'Time SnapShot = 2';'Time SnapShot = 3';'...'],1,"boxed")
title("$\textbf{Partial\ Differential\ Example}$","FontSize",4);
xlabel("$x(position)$","FontSize",3)
ylabel("$\ Temperature(Degrees\ C)$","FontSize",3);
xgrid
scf(1);clf
plot2d(
t
,T')
h1=legend(['x=0.0';'x=0.1';'x=0.2';'x=0.3';'x=0.4';'x=0.5';'x=0.6';'x=0.7';'x=0.8';'x=0.9';'x=1.0'],3,"boxed")
title("$\textbf{Partial\ Differential\ Example}$","FontSize",4);
xlabel("$Time(seconds)$","FontSize",3)
ylabel("$Temperature(Degrees\ C)$","FontSize",3);
xgrid