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