MATLAB Answers

Solving Finite Difference Method Using ODE15s

28 views (last 30 days)
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.

  0 Comments

Sign in to comment.

Accepted Answer

Torsten
Torsten on 8 Aug 2019
Edited: Torsten on 8 Aug 2019
%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).

  3 Comments

Mathushan Sabanayagam
Mathushan Sabanayagam on 8 Aug 2019
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
Torsten on 9 Aug 2019
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
Mathushan Sabanayagam on 9 Aug 2019
Thank you for explaining, it makes sense. Thank you for all your help, really appreciate it!

Sign in to comment.

More Answers (1)

Torsten
Torsten on 7 Aug 2019
Edited: Torsten on 7 Aug 2019
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.

  1 Comment

Mathushan Sabanayagam
Mathushan Sabanayagam on 7 Aug 2019
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)

Sign in to comment.

Sign in to answer this question.

Tags