Problem to plot the results of a differential equation solved with ode15s

Hello everyone,

I need help with my code. I want to plot the degree of cure of epoxy as a function of time for several temperature using the Khoun model.

First, I have to solve a differential equation. To do that, I use ode15s. Then, I want to plot my results.

My problem is that matlab tells me that there is no error in my code. However, The results displayed are wrong, the curves displayed are horizontal lines starting at zero.

My code :

% This script calculate values of the degree of cure based on the Khoun model
clear
close all

% Fitting parameters values for the cure kinetics

A = 481955; % s^-1
E = 54591; % Jmol^-1
n = 1.87;
m = 0.05;
C = 61.7;
a_c = 0.97;
a_T = 0.0014

; % K^-1
R = 8.314;

tmax = 500;
n_point = 10;
tspan = linspace(0,tmax,n_point);

a0 = 0; % Initial condition of alpha

T_ref = [50 70 90 110]; % Temperature, °C

for iter = 1:1:length(T_ref)
    T = T_ref(iter);

    [t,a] = ode15s(@(t,a) ((A*exp(-E/(R*(T+273.15))))/(1+exp(C*(a-a_c-a_T*(T+273.15)))))*((1-a)^n)*a^m, tspan, a0); % Solving the equation of the degree of cure
    
    for i = 1:1:n_point
        a = abs(a);    % The imaginary part can be neglect
        if a(i) > 1 
        a(i) = 1;   % It's not possible to have a degree of cure higher than 1
        end
    end
    
    hold on;
    plot(t,a);
    xlabel('Time, (min)');
    ylabel('Degree of cure \alpha, (-)');
    title('Evolution of the degree of cure \alpha versus time - Khoun cure kinetic model');
    legend('50°C','70°C','90°C','110°C','Location','Best');
end

Here is what is plotted which is wrong.

Do you have an idea of what is wrong with my code ?

Thank you for your help.

Faouzi

Answer 1

[Star Strider]

Use a value for your initial condition that is near 0, not equal to 0:

a0 = 1E-15; % Initial condition of alpha

That produces the family of curves you likely want.

[Faouzi Gheffar]

Thank you for your help :)

[Star Strider]

As always, my pleasure!