Help pls, i have problem with using for loop to plot 2 input variables equation.

A -k1-> B -k2-> C

A -k3-> D

F and Qdot are input variables.

clear all
%Parameter
K0ab = 1.287*(10^12);        %h^-1
K0bc = 1.287*(10^12);        %h^-1
K0ad = 9.043*(10^9);         %l/mol*h
Rgas = 8.3144621*(10^(-3));
EAab = 9758.3;               %kJ/mol
EAbc = 9758.3;               %kJ/mol
EAad = 8560.0;               %kJ/mol
HRab = 4.2;                  %kJ/molA
HRbc = -11.0;                %kJ/molB
HRad = -41.85;               %kJ/molA
Tin = 130.0;                 %degree C
Kw = 4032.0;                 %kJ/(h*(m^2)*K)
Rou = 0.9342;                %kg/l
Cp = 3.01;                   %kJ/kg*K
Cpk = 2.0;                   %kJ/kg*K
AR =0.215;                   %m^2
VR = 10.01;                  %l
mk = 5.0;                    %kg
%Initial condition
CA(1) = 5.1;                 %mol/l = CA0
CB(1) = 0.8;                 
TR(1) = 140;                 
TK(1) = 140;                 
%Step size
dt=0.01;
tt=50;
n = tt/dt ;
t(1) = 0 ;

%System simulation
for i = 1:n
    for F = 5:1:100          %between(5,100)
        for Qdot = -8500:1:0    %between(-8500,0)
            k1 = K0ab*exp((-EAab)/(TR(i)+273.15));
            k2 = K0bc*exp((-EAbc)/(TR(i)+273.15));
            k3 = K0ad*exp((-EAad)/(TR(i)+273.15));
            CA(i+1) = CA(i)+dt*((F*(CA(1)-CA(i)))-(k1*CA(i))-(k3*(CA(i)^2)));
            CB(i+1) = CB(i)+dt*((-(F*CB(i)))+(k1*CA(i))-(k2*CB(i)));
            TR(i+1) = TR(i)+dt*((((k1*CA(i)*HRab)+(k2*CB(i)*HRbc)+(k3*(CA(i)^2)*HRad))/(-(Rou*Cp)))+(((F*(Tin-TR(i)))+(Kw*AR*(TK(i)-TR(i))))/(Rou*Cp*VR)));
            TK(i+1) = TK(i)+dt*((Qdot+(Kw*AR*(TK(i)-TR(i))))/(mk*Cpk));
            t(i+1) = t(i) + dt;
        end
    end
end
F = 5:1:100;          %between(5,100)
Qdot = -8500:1:0;     %between(-8500,0)

%Plot graph
figure
subplot(3,2,1)
plot(t,CA,'-')
xlabel('time')
ylabel('C_A')
subplot(3,2,2)
plot(t,CB,'-')
xlabel('time')
ylabel('C_B')
subplot(3,2,3)
plot(t,TR,'-')
xlabel('time')
ylabel('T_R')
subplot(3,2,4)
plot(t,TK,'-')
xlabel('time')
ylabel('T_K')
subplot(3,2,5)
plot(t,F,'-')
xlabel('time')
ylabel('F')
subplot(3,2,6)
plot(t,Qdot,'-')
xlabel('time')
ylabel('Qdot')

[Torsten]

Use SI units. At the moment, you work with h and J which is not compatible.

Further, you vary three parameters, namely time, F and Qdot, Therefore, CA, CB, TR and TK must be three-dimensional arrays: CA(time,F,Qdot), CB(time,F,Qdot), TR(time,F,Qdot) and TK(time,F,Qdot).

Answer 1

[Ruchika Parag]

Hi @athittaya, you're facing issues because you're looping over F and Qdot inside the simulation loop without storing results separately. This causes overwriting and incorrect plots. Also, F and Qdot need to be fixed values during integration unless you're running parameter sweeps.

Here's a working example that simulates the system for a single pair of inputs: F = 50, Qdot = -5000, with dummy values added for missing parameters:

clear all
% Constants and dummy parameters
K0ab = 1.287e12; K0bc = 1.287e12; K0ad = 9.043e9;
EAab = 9758.3; EAbc = 9758.3; EAad = 8560.0;
HRab = 4.2; HRbc = -11.0; HRad = -41.85;
Tin = 130.0; Kw = 4032.0; Rou = 0.9342;
Cp = 3.01; Cpk = 2.0; AR = 0.215; VR = 10.01; mk = 5.0;
% Time setup
dt = 0.01; tt = 50; n = tt/dt;
% Inputs
F = 50; Qdot = -5000;
% Initial conditions
CA = zeros(1,n+1); CB = CA; TR = CA; TK = CA; t = CA;
CA(1) = 5.1; CB(1) = 0.8; TR(1) = 140; TK(1) = 140;
% Simulation loop
for i = 1:n
    k1 = K0ab * exp(-EAab / (TR(i)+273.15));
    k2 = K0bc * exp(-EAbc / (TR(i)+273.15));
    k3 = K0ad * exp(-EAad / (TR(i)+273.15));
    CA(i+1) = CA(i) + dt * ((F*(CA(1)-CA(i))) - k1*CA(i) - k3*(CA(i)^2));
    CB(i+1) = CB(i) + dt * ((-F*CB(i)) + k1*CA(i) - k2*CB(i));
    TR(i+1) = TR(i) + dt * ((((k1*CA(i)*HRab) + (k2*CB(i)*HRbc) + (k3*CA(i)^2*HRad))/(-(Rou*Cp))) ...
               + ((F*(Tin - TR(i)) + Kw*AR*(TK(i) - TR(i))) / (Rou*Cp*VR)));
    TK(i+1) = TK(i) + dt * ((Qdot + Kw*AR*(TK(i)-TR(i))) / (mk*Cpk));
    t(i+1) = t(i) + dt;
end
% Plot
figure
subplot(2,2,1), plot(t, CA), xlabel('Time'), ylabel('C_A')
subplot(2,2,2), plot(t, CB), xlabel('Time'), ylabel('C_B')
subplot(2,2,3), plot(t, TR), xlabel('Time'), ylabel('T_R')
subplot(2,2,4), plot(t, TK), xlabel('Time'), ylabel('T_K')

●

Once this works, you can extend it for a sweep across F and Qdot using nested loops and storing results accordingly. Thanks!