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

Question

Faouzi Gheffar 2019-4-12

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/456024-problem-to-plot-the-results-of-a-differential-equation-solved-with-ode15s

评论： Star Strider 2019-4-12

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