MATLAB Answers

ODE for the modelling of polymers degradation - products of sums, additional matrix in the function?

조회 수: 4(최근 30일)
Aga
Aga 2021년 5월 12일
댓글: Bjorn Gustavsson 2021년 6월 17일
EDIT TO ADD: The below approach to depolymerization is wrong!
Hello MATLAB central,
I am trying to model a sequence of bond fission and radical recombination reactions in a polymer. In the below figure reaction (1) is the bond fission reaction and reaction (2) is the radical recombination reaction.
(Figure adapted from Levine et.al, 2009)
I am using kinetic equations of the form:
(1) where n is the length of the polymer A (of concentration [An]), [Rn] are the concentration of radicals R of length n, k1 and k2 are the reaction rates of the processes
(2)
I am planning to use an ODE solver to get the concentration of all evolved species, but I am struggling with defining the sum of the products in reaction (1).
I first wanted to introduce additional matrix S that would describe products of all possible combinations of radicals that could produce the alkane of length j (see my code below), but it doesn't work. When using a stiff solver It returns a matrix of M with a warning:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.004156e-20.
While non-stiff solvers produce no warning message. The resulting matrix is just the initial conditions row repated for length of the timespan.
The code is: (please note that I used 23 as the length of the initial polymer, and I know I should probably have defined it as an additional argument of this function and I will change it after figuring out this issue).
% Levine model
function dMdt = levine(M,~,~, k1, k2, sumS, S)
dMdt = zeros(45,1);
for n = 1:45
for i =1:22 %additional array to account that not all radicals can form all alkanes
if n>=24 && n<=45 && i==1 % range of n's that are corresponding to the location of radicals and special case (i=1) for
% the formation of alkane of length 23
S(i,n) = M(n).*M(69-n); %all of the possible combination of radicals that can form alkane of length 23
elseif n >=24 && n<=46-i && i>1 && i<=22 %range of the n's that correspond to the location of radicals in the array that can
% combine to create alkanes of length 2 to 22
S(i,n) = M(n).*M(70-i-n); %all of thecombination of radicals that can combine to create alkane of length 2:22
sumS = sum(S);
if n>=1 && n<=22
dMdt(n) = -k1*M(n)+0.5*k2*sumS(n); %alkanes (excluding methane);
elseif n==23 %methane
dMdt(n) = 0; %does not undergo fission and isnt generated through radical recombination
elseif n>23 && n<=45 %range of primary radicals in the array
dMdt(n) = -k1*M(n)+k2*sum(M(n-22:24));
end
end
end
end
And the command that I use to call this function is:
tspan = 1:1800;
m0 = [1 zeros(1,44)];
T = 673;
k1 = 1*10^16*exp(-89.7/(T*1.987*10^-3)) ;
k2 = 1.10*10^11*exp(-2.3/(T*1.987*10^-3)) ;
[t,M] = ode23s(@(t,M) levine(M,n,i,k1,k2,S,sumS), tspan, m0);
  댓글 수: 3
Bjorn Gustavsson
Bjorn Gustavsson 2021년 6월 17일
Ah, that long molecules and that messy mixtures. I can see that that might become challenging. The chemistry I manage only contain some 10 reacting species (O, O2, N, N2 and a few of their ions and molecules) and there we can brute-force it.

댓글을 달려면 로그인하십시오.

채택된 답변

Bjorn Gustavsson
Bjorn Gustavsson 2021년 5월 12일
It seems to me that you might have your problems when calculating the products of the radicals. Your ODE-function for the coupled odes of the radicals and alkenes might look something like this:
function dARdt = polymer_reactions(t,AR,k1,k2)
A = AR(1:end/2);
R = AR((1+end/2):end); % simple partitioning here, because I'm not a chemist
dARdt = [-k1*A(:)
-k2*R(:)];
for n = 1:numel(A) % explicitly loop over all combinations
for j1 = 1:(n-1) % not necessarily elegant as such...
dARdt(n) = dARdt(n) + k2*R(j1)*R(n-j1);
end
end
for n = 1:numel(R)
dARdt(end/2+n) = dARdt(end/2+n) + k1*sum(A((n+1):end)); % This seems a bit idealized, no?
end % not different rates to produce radicals of different lengths from alkenes of any length?
HTH
  댓글 수: 3

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by