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.

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).

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))

Sign in to comment.

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.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.