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

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

[Bjorn Gustavsson]

In what sense is it "wrong"? Not correctly taking into how the polymers partition into different smaller chains?

[Aga]

Hello Bjorn,

After having researched the issue further, I realized that the above approach would work if I wanted to monitor each of the species that are possible to be evolved from the polymer decomposition, which would make the model have to include thousends of reactions! (I misunderstood the initial paper by a mile :( ). Also, the initial polymer is not one long chain, but rather a mixture of chains of different lengths. The polymers are usually represented by their molecular weight distribution.

Currently, I am working on a statistical method that predicts the changes of moments of that distribution during pyrolysis.

[Bjorn Gustavsson]

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.

Answer 1

[Bjorn Gustavsson]

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

[Aga]

Thank you very much Bjorn!

Re: rates of production of radicals of different sizes from alkanes: this model assumes the same reaction rate for all reactions of the same type, so thankfully I don't have to worry about that :).

[Bjorn Gustavsson]

Good that it helped.

About chemistry: Figures, but simplifying models are either good enough, or a (nice) starting-point, right?

[Aga]

Very true!