Unable to perform assignment because the left and right sides have a different number of elements.
1 次查看(过去 30 天)
显示 更早的评论
Hi all, I don't know why does my code give an error message:
Unable to perform assignment because the left and right sides have a different number of elements.
Error in eulerian>odesfun (line 115)
dydt(1) = (1 ./V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2 ;
Error in eulerian (line 53)
dydt = odesfun(t,V_Headspace, F, CSO2_in, CCO2_in, NSO2, NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g,y0);
The code is given below, and on the attachment:
function eulerian()
global S_total C_total Ca_total CSO2_bulkslurry CSO3 ...
CCO2_bulkslurry CCaCO3 CCaSO3 CH CCO2_in...
E NSO2 E_CO2 NCO2 RCaSO3 RCaCO3 SO2_g CO2_g
global V_Headspace F CSO2_in R Temp HSO2 kga kLa DCa2 ...
DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 ...
MWCaCO3 Kad KCO2 KHCO3 KSO2 KHSO3 Kw CH0
% Constant Parameters
KSO2 = 6.24;
KHSO3 = 5.68e-5;
DCa2 = 1.39e-9;
DSO2 = 2.89e-9;
R = 8.314;
Temp = 323.15;
HSO2 = 149;
kga = -4.14e-6;
kLa = -7e-2;
V_Headspace = 1.5e-6;
F = 1.67e-4;
CSO2_in = 6.51e-2;
KCO2 = 1.7e-3;
KHCO3 = 6.55e-8;
DCO2 = 3.53e-9;
kLa_CO2 = 9.598e-4;
HCO2 = 5.15e+3;
CCO2_in = 0;
KSPCaSO3 = 1.07e-7;
BETCaSO3 = 10;
ktot = 8.825e-3;
BETCaCO3 = 12.54;
MWCaCO3 = 100.0869;
Kad = 0.84;
CH0 = 7.414e-6;
Kw = 5.3e-8;
tspan = [600;1200;1800;2400;3000;10200;17400;24600;31800;39000;...
46200;53400;60600;67800;75000;82200;89400;96600;103800;111000;...
118200;125400;132600;139800;147000;154200;161400;168600;175800;183000];
y0 = [1.72E-02; 6.97E-08; 3.34E-01; 4.88E-06; 4.88E-06; 4.86E+01; 0];
y0 = y0';
delta = 0.01;
nloop = 30;
t0 = 600.;
t = t0;
for i = 1:nloop
dydt = odesfun(t,V_Headspace, F, CSO2_in, CCO2_in, NSO2, NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g,y0);
y = y0+ dydt.'*delta;
%Variables SO2_g, CO2_g, S_total, C_total, Ca_total, CCaCO3 and CCaSO3 are passed as y(1)...y(7)
SO2_g = y(1);
CO2_g = y(2);
S_total = y(3);
C_total = y(4);
Ca_total = y(5);
CCaCO3 = y(6);
CCaSO3 = y(7);
% CH, CSO2_bulkslurry , E, NSO2, CCO2_bulkslurry, E_CO2, NCO2, CSO3, RCaSO3 and RCaCO3 are passed as variable Parameters
CH = fsolve(@(CH) CH + 2.* Ca_total - ((S_total.* KSO2.*CH)/(CH^2 + KSO2.*CH + KSO2.*KHSO3))- 2.*((S_total.*KSO2.*KHSO3)/(CH^2 ...
+ KSO2.*CH + KSO2.*KHSO3))-((C_total.*KCO2.*CH)/(CH^2 + KCO2.*CH + KCO2.*KHCO3))- 2.*((C_total.*KCO2.*KHCO3)/(CH^2 ...
+ KCO2.*CH + KCO2.*KHCO3))- Kw/CH, CH0);
CSO2_bulkslurry = (S_total.* CH.^2)/(CH.^2 + KSO2.* CH + KSO2.* KHSO3);
E = 1 + DCa2.* Ca_total / (DSO2.* S_total);
NSO2 =(SO2_g.* R.* Temp - HSO2.* CSO2_bulkslurry) / ( (1/kga) + HSO2/ (E.* kLa) ) ;
CCO2_bulkslurry = (C_total.* CH.^2)/(CH.^2 + KCO2.* CH + KCO2.* KHCO3);
E_CO2 = 1 + DCa2 * Ca_total / (DCO2 .* C_total);
NCO2 = kLa_CO2 .* E_CO2 .* (CO2_g .* R .* Temp / HCO2 - CCO2_bulkslurry);
CSO3 = (S_total.* KSO2 .* KHSO3)/(CH.^2 + KSO2 .* CH + KSO2 .* KHSO3);
RCaSO3 = 0.001 .* 1.62e-3 .* exp(-5152/ Temp) .* ( (Ca_total .* CSO3 / KSPCaSO3) - 1).^3 / BETCaSO3;
RCaCO3 = - ktot .* BETCaCO3 .* MWCaCO3 .* CCaCO3 .* CH .* (1 - (Kad .* CH)/(1+ Kad .* CH));
y0 = y;
t = t + delta;
CH0 = CH;
[t,y] = ode45(@odesfun,tspan,y0)
figure (1)
plot (t,y)
legend('Mod-SO_{2-gas}', 'Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}',...
'Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^+','Mod-pH', 'Location','best')
end
end
function dydt = odesfun(t,y,V_Headspace, F, CSO2_in, CCO2_in, NSO2, NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g)
%Variables SO2_g, CO2_g, S_total, C_total, Ca_total, CCaCO3 and CCaSO3 are passed as y(1)...y(7)
dydt = zeros(7,1);
% SO2 in gas phase
dydt(1) = (1 ./V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2 ;
% CO2 in gas phase.
dydt(2) = (1 ./V_Headspace).* (F .* CCO2_in - F .* CO2_g) - NCO2 ;
% Total sulfur balance
dydt(3) = NSO2 - RCaSO3;
% Total sulfur balance
dydt(4) = RCaCO3 - NCO2;
% Total calcium balance
dydt(5) = RCaCO3 - RCaSO3;
% Limestone dissolution submodel
dydt(6) = RCaCO3;
% Hannebachite crystalization submodel
dydt(7) = RCaSO3;
end
Please help.
7 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Oceanography and Hydrology 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!