Solving Finite Difference Method Using ODE15s

Hi,

I have created a function in which I am solving a partial differential equation where temperature is dependent on time and radius (energy balance in spherical coordinates). I have discretised the spatial coordinates into n nodes. This yielded an ode with respect to time that I must complete over each node. There are also two other differential equations (not relevant to this however, its solutions are stored in n+1 and n+2 of the DyDt matrix just for information). I have shown the relevant parts of the code for brevity:

T = zeros(n,1); %initialise T as a matrix

DTDt = zeros(n,1); %initialise DTdt as a matrix

DyDt = zeros(n+2,1); %initialise a matrix containing DTdt from 1:n. (the n+1 and n+2 are the two other solutions)

T = y(1:n); %fill the T values into a y matrix from 1 to nth column

I then have a for loop for i=1:n-1 with my expression for DTDt. This allows me to solve all nodes up unti n-1. However, my problem is for the solution to T(n), I do not have a DTdt equation, but instead I have an algebraic equation:

T(n)=((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr));

It relies only only values from the other nodal positions and some other constants.

In a separate file, I utilise the functions, setting the initial condition as T0 at 298K.

How can I go about solving such a system? I understand it is a set of differential equations and 1 algebraic equation.

Thank you very much in advance. I am a beginner in Matlab and will appreciate any help.

Answer 1

[Torsten]

%OTHER FILE THAT CALLS FUNCTION (Excluding all the constants)
T0 = ones(n,1)*298; %matrix of initial condition for T, rows = number of nodes, only 1 column since t=0
a0 = 1; %initial condition for a
b0 = 1; %initial condition for b
y0 = [T0;a0;b0]; %T0 followed by a0 and b0 in a single column matrix 
M = eye(n+2);
M(n,n) = 0.0;
options = odeset('Mass',M)
%Call the solver
[T Y] = ode15s(@(t,y)convection(t,y,n,preexpomelt,activationmelt,preexpodecom,enthalpymelt,enthalpydecom,activationdecom,heatrate,alphas,alphal,dr,cps,R0,cpl,Tm,h,Tg,ks),t,y0,options)

And in "convection", set

dTdt(n) = T(n)-(((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr)));

This is the second method which should be easier to implement since the solver still has the same number of unknowns as before (n+2).

[Mathushan Sabanayagam]

Hi Torsten,

Thanks for the help! I was wondering why the T(n)=((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr)) equation is now replaced by dTdt(n) = T(n)-(((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr))); ? Also, the initial condition of T(n)=0 does not seem to work anymore

Thank you for all your help!

[Torsten]

I was wondering why the T(n)=((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr)) equation is now replaced by dTdt(n) = T(n)-(((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr)));

Since the element in the mass matrix corresponding to dTdt(n) is set to zero, the equation now reads

0*dTdt(n) = T(n)-(((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr))); 

Also, the initial condition of T(n)=0 does not seem to work anymore

The initial condition T0 for T(n) should be

T0 = ((2.*h.*Tg.*dr)-(ks.*T0(n-2))+(4.*ks.*T0(n-1)))/((3.*ks)+(2.*h.*dr))

[Mathushan Sabanayagam]

Thank you for explaining, it makes sense. Thank you for all your help, really appreciate it!

Answer 2

[Torsten]

You can either solve for T1,...,T_(n-1) (and the two extra variables) using ODE15S and calculate

T(n) internally each time the function to calculate the time derivatives is called or you can use the mass matrix option of ODE15S by setting M(n,n) = 0 instead of M(n,n) = 1. This will tell ODE15S that equation n is an algebraic equation instead of a differential equation.

[Mathushan Sabanayagam]

Hi Torsten,

Thank you for your reply! Sorry if its an obvious question but how would I go about calculating the T(n) internally. Would the expression be in the function file or the file that calls the function? I would prefer the first method as I am already using ODE15s and so I would find it easier to implement. Thanks in advance for your help!

I've attached the code if that helps:

function DyDt=convection(t,y,n,preexpomelt,activationmelt,preexpodecom,enthalpymelt,enthalpydecom,activationdecom,heatrate,alphas,alphal,dr,cps,R0,cpl,Tm,h,Tg,ks)

%INITIALISING VARIABLES AND DERIVATIVES AS MATRICES
T = zeros(n,1); 
a=0; 
b=0;
DTDt = zeros(n,1); 
DaDt = 0; 
DbDt=0; 
DyDt = zeros(n+2,1);
DTDr = zeros(n-1,1); 
D2TDr2 = zeros(n-1,1); 

%STORING SOLUTIONS IN Y MATRIX
T = y(1:n); %fill the T values into a y matrix from 1 to nth element
a = y(n+1); %fill the a values into a y matrix at the n+1th element 
b= y(n+2);  %fill the a values into a y matrix at the n+2th element 

%TEMPERATURE-RADIUS DERIVATIVE     
for i=2:n-1 
  DTDr(i) = ( T(i+1) - T(i) )./ (dr) ; 
  D2TDr2(i) = ( 1./(dr^2) ) .* (  T(i+1) - (2.*T(i)) + T(i-1) ); 
end

%TEMPERATURE TIME DERIVATIVE 
for i=2:n-1 
  DTDt(i) = ((((((a.^3)/(b.^3)).*alphas) +  ((((b.^3)-(a.^3))/b.^3).*alphal))./(((R0.*b).^2))).* (D2TDr2(i))) + ((((((a.^3)./(b.^3)).*alphas) +  ((((b.^3)-(a.^3))./b.^3).*alphal))./(((R0.*b).^2))).* (DTDr(i).*(2./(i.*dr))))...
      +  (((1./((((a.^3)./(b.^3)).*cps) +((((b.^3)-(a.^3))./b.^3).*cpl)))).*((enthalpymelt.*preexpomelt.*exp(-activationmelt./(8.314.*T(i))))+(enthalpydecom.*preexpodecom.*exp(-activationdecom./(8.314.*T(i))))));
end
 
for i=1
    DTDt(i)=DTDt(i+1);
end

T(n)=((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr));

if T(n) < Tm 
DaDt(1)=0;
else 
DaDt(1)= ((-1.*a.*preexpomelt)./3).*exp(-activationmelt./((8.314.*Tm)));
end
DbDt(1)= ((-1.*b.*preexpodecom)./3).*exp(-activationdecom./((8.314.*T(n))));  
DyDt = [DTDt;DaDt;DbDt]; 
end

%OTHER FILE THAT CALLS FUNCTION (Excluding all the constants)
T0 = ones(n-1,1)*298; %matrix of initial condition for T, rows = number of nodes, only 1 column since t=0
a0 = 1; %initial condition for a
b0 = 1; %initial condition for b
y0 = [T0;a0;b0]; %T0 followed by a0 and b0 in a single column matrix 
%Call the solver
[T Y] = ode15s(@(t,y)convection(t,y,n,preexpomelt,activationmelt,preexpodecom,enthalpymelt,enthalpydecom,activationdecom,heatrate,alphas,alphal,dr,cps,R0,cpl,Tm,h,Tg,ks),t,y0)