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

2 次查看（过去 30 天）
athittaya
athittaya 2025-4-29
回答： Ruchika Parag 2025-6-5
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')
  1 个评论
Torsten
Torsten 2025-4-29
编辑：Torsten 2025-4-29
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).

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Ruchika Parag
Ruchika Parag 2025-6-5
Ran in:
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!

类别

在 Help Center 和 File Exchange 中查找有关 Genomics and Next Generation Sequencing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by