Solving Finite Difference Method Using ODE15s

15 次查看（过去 30 天）
Mathushan Sabanayagam
评论： Mathushan Sabanayagam 2019-8-9
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.

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2019-8-8
编辑：Torsten 2019-8-8
%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 个评论
显示 1更早的评论
Torsten
Torsten 2019-8-9
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 2019-8-9
Thank you for explaining, it makes sense. Thank you for all your help, really appreciate it!

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2019-8-7
编辑：Torsten 2019-8-7
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 个评论
Mathushan Sabanayagam
Mathushan Sabanayagam 2019-8-7
编辑：Mathushan Sabanayagam 2019-8-7
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)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by