MATLAB Answers

How to use fprintf

Dursman Mchabe 님이 질문을 제출함. 22 Jan 2019
최근 활동 per isakson 님이 편집함. 22 Jan 2019
How can I use fprintf to print out
NSO2;
CSO2_bulkslurry;
CHSO3;
CSO3;
CSO2_inter;
CHSO3_inter;
CSO3_inter
on the script below:
global S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 Trial CH P CCO2_in nloop Delt t CH_trial ...
CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter CCO3_inter pHtrial pi pHInter V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad ...
rhoCaCO3 KCO2 KHCO3 KSO2 KHSO3 Kw MWCaSO3 rhoCaSO3 kCaS03 PercentageAbsorption DHSO3 DSO3 Enhancement NSO2
format shorte
%% Parameters
%Symbol Value Description Units
R = 8.314; % Universal gas constant (J/molK)
HSO2 = 149; % SO2 Henry’s law constant (m3Pa/mol)
HCO2 = 5.15e+3; % CO2 Henry’s law constant (m3Pa/mol)
DCa2 = 1.39e-9; % Diffusivity of Ca2+ ions (m2/s)
DSO2 = 2.89e-9; % Diffusivity of SO2 (m2/s)
DCO2 = 3.53e-9; % Diffusivity of CO2 (m2/s)
KCO2 = 1.7e-3; % CO2 dissociation equilibrium constant (mol/m3)
KHCO3 = 6.55e-8; % HCO3- dissociation equilibrium constant(mol/m3)
KSO2 = 6.24; % SO2 dissociation equilibrium constant (mol/m3)
KHSO3 = 5.68e-5; % HSO3- dissociation equilibrium constant(mol/m3)
Kw = 5.3E-8 ; % water dissociation constant [mol/m3.]^2
KSPCaSO3 = 3.1e-4; % Solubility product of CaSO3 (mol2/m6)
BETCaSO3 = 10; % BET specific surface area of CaSO3 (m2/g)
ktot = 8.825e-3; % Total dissolution rate constant (m/s)
BETCaCO3 = 12.54; % BET specific surface area of CaSO3 (m2/g)
MWCaCO3 = 100.0869; % Molecular weight of CaCO3 (g/mol)
Kad = 8.4e-3; % Adsorption constant (m3/mol)
rhoCaCO3 = 2703; % Density of limestone (kg/m3)
MWCaSO3 = 258.30; % Molecular weight of CaSO3.0.5H2O (g/mol)
rhoCaSO3 = 2540; % Density of CaSO3.0.5H2O (kg/m3)
kCaS03 = 2e-08; % rate constant for CaSO3 (mol/m^3 s)
P = 1.01E+05; % operating pressure (Pa)
T = 323.15; % Reactor temperature (K)
V_Headspace = 4.e-4; % Volume of the headspace (m3)% changed t0 -4
F = 1.66667e-5; % Flue gas inlet flow rate (m3/s)
y_so2 = 2000* 1E-6; % 2000 ppm
%CSO2_in = y_so2 * P/R/T;
CSO2_in = 7.28E-02; % Inlet flue gas SO2 concentration (mol/m3)
y_co2 = 0.0;
CCO2_in = y_co2 * P/R/T;% CO2 concentration in gas phase, (molm^3)
kga = 4.14e-6; % Product of gas-side mass coefficient and interfacial surface area (mol/m3/Pa/s)
kga = 4.74e-5; % suggested value
kLa = 6e-4; % Product of liquid-side mass coefficient and interfacial surface area for SO2 (1/s)
kLa_CO2 = 9.598e-5; % Product of gas-side mass coefficient and interfacial surface area for SO2 (1/s)
DHSO3 = 2.04E-09; % Diffusivity of HSO3 (m2/s)
DSO3 = 1.70E-09; % Diffusivity of SO3 (m2/s)
%% Initialization
pHtrial = 6.8;
Ca_total = 0 * 4.879E-02;
C_total = 0 * 4.879E-02;
S_total = 0.0;
CCaCO3 = 0.0;
CCaSO3 = 0.0;
% SO2_g is a fraction (fr) of CSO2_in
fr = 0.01;
SO2_g = fr * CSO2_in;
%
CO2_g = CCO2_in;
CO2_g = 0.2;
CCO2_inter = CO2_g ;
%
pg = 200;
CSO2_inter = pg/HSO2;
% Initial pH (bulk and interface)
CH_trial = 1.e-7;
pH1 = 6.8;
pH = fzero (@HionpH, pH1 );
CH = 10^(3-pH) ;
%% C_Bulk
% Relationship between CSO2_bulkslurry, S_total,CH+,KSO2 and KHSO3-
CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between HSO3-, S_total,CH+,KSO2 and KHSO3-
CHSO3 = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between SO32-, S_total,CH+,KSO2 and KHSO3-
CSO3 = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between CCO2_bulkslurry, C_total,CH+,KCO2 and KHCO3-
CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);
% Relationship between HCO3-, C_total,CH+,KCO2 and KHCO3-
CHCO3 = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);
% Relationship between CO32-, C_total,CH+,KCO2 and KHCO3-
CCO3 = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
%%
pH2 = 6.0;
pHInter = fzero (@HionStar, pH2);
CH = 10^(3-pH) ;
% set initial conditions
y0 = [ SO2_g CO2_g S_total C_total Ca_total CCaCO3 CCaSO3 ];
y0 = y0'; % put as a column vector
% set trial values for Fsolve eo be used in ODES
Trial(1) = (SO2_g * R *T ) / HSO2; % assumed so2 at interface
Trial(2) = 1; % SO2_total at interface
Trial(3) = CH ; % hydrogen at interfafe
%% time stepping parametes
Delt = 1.0 ; % time step changed from 0.01;
t0 = 0; % time for further simulation
t = t0;
tmax = 930.90909090909090909090;% maximum time for simulation for each loop.
nloop = tmax/Delt ; % number of steps to be taken.
nstep = 33 ;
% template for writing the solution.
Results = [ t y0' pH pHInter];
%% time stepping starts.
for k = 1:nstep
% loop for nloop time steps.
y = SO2_OdeDriver( y0);
% write the results and repeat.
pH = 3 - log10 (CH) ;
tf = t;
% intermediate results
Results = [Results; tf y' pH pHInter];
PercentageAbsorption = ((CSO2_in - Results(:,2))/CSO2_in).*100;
Enhancement = 1 + DHSO3.* (CHSO3_inter - CHSO3)/DSO2.* (CSO2_inter - Results(:,2)) + DSO3.* (CSO3_inter - CSO3)/ DSO2.* (CSO2_inter - Results(:,2));
%NSO2 = (SO2_g* R* T - HSO2* Results(:,2))./((1/kga) + HSO2/(Enhancement* kLa));
y0 = y;
end
Results;
Solution = [ Results(:,1) Results(:,2) Results(:,4) Results(:,9)];
PercentageAbsorption = ((CSO2_in - Results(:,2))/CSO2_in).*100
Enhancement
Exp_Time=[0;960;1920;2880;3840;4800;5760;6720;7680;8640;9600;10560;11520;12480;13440;14400;15360;16320;17280;18240;19200;20160;21120;22080;23040;24000;24960;25920;26880;27840;28800;29760;30720];
Exp_SO2_g =[7.28E-02;3.66E-02;4.13E-02;4.44E-02;4.78E-02;5.02E-02;5.30E-02;5.48E-02;5.70E-02;5.88E-02;6.03E-02;6.18E-02;6.35E-02;6.43E-02;6.52E-02;6.62E-02;6.71E-02;6.81E-02;6.85E-02;6.92E-02;6.98E-02;7.02E-02;7.06E-02;7.11E-02;7.12E-02;7.15E-02;7.18E-02;7.20E-02;7.22E-02;7.23E-02;7.25E-02;7.26E-02;7.27E-02];
Exp_S_total = [0;0.35;0.69;0.99;1.26;1.51;1.72;1.91;2.08;2.23;2.37;2.49;2.59;2.68;2.76;2.83;2.89;2.94;2.99;3.03;3.06;3.09;3.11;3.13;3.15;3.16;3.17;3.18;3.19;3.19;3.20;3.20;3.20];
Exp_pH = [6.87;6.65;6.43;5.68;3.45;3;2.85;2.78;2.63;2.61;2.57;2.56;2.55;2.58;2.59;2.55;2.57;2.57;2.52;2.51;2.51;2.53;2.53;2.51;2.5;2.51;2.49;2.49;2.49;2.47;2.47;2.46;2.46];
figure (1)
plot(Exp_Time,Exp_SO2_g,'kd',Results(:,1), Results(:,2),'k-')
legend('Exp-SO_{2,g}','Mod-SO_{2,g}')
ylabel('Conc (mol/m^{3}','FontSize',12,'FontWeight','bold')
xlabel ('Time (sec)','FontSize',12,'FontWeight','bold')
set(gca,'YColor','k');
figure(2)
plot(Exp_Time,Exp_S_total,'bo', Results(:,1),Results(:,4),'b-')
legend('Exp-S_{total}','Mod-S_{total}','location','best')
ylabel('Conc (mol/m^{3}','FontSize',12,'FontWeight','bold')
xlabel ('Time (sec)','FontSize',12,'FontWeight','bold')
set(gca,'YColor','k');
figure(3)
plot(Exp_Time,Exp_pH,'rv',Results(:,1),Results(:,9),'r-')
legend('Exp-pH','Mod-pH','location','best')
ylabel('pH ','FontSize',12,'FontWeight','bold')
xlabel ('Time (sec)','FontSize',12,'FontWeight','bold')
set(gca,'YColor','k');
%% pHmodel
function ph = HionpH (pH)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 ...
KSO2 KHSO3 Kw MWCaSO3 rhoCaSO3 S_total C_total Ca_total
% Ph suppied. calculate CH and set up discrepacncy function.
CH = 10^(3-pH);
% Hydrogen balance from electro-neutrality.
ph = 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 ;
end
function ph1 = HionStar (pH)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 ...
KSO2 KHSO3 Kw MWCaSO3 rhoCaSO3 S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter ...
CCO3_inter CCO2_bulkslurry CHCO3 CCO3 pg
CH = 10^(3-pH);
% species concentrations.
CSO2_inter ;
CHSO3_inter = KSO2*CSO2_inter / CH ;
CSO3_inter = KHSO3 *CHSO3_inter /CH;
Stotal = CHSO3 + CSO3;
Ctotal = CHCO3 + CCO3 ;
CHSO3_inter = Stotal * ( 1+ KHSO3 / CH) ;
CSO3_inter = KHSO3 * CHSO3_inter /CH;
CHCO3_inter = Ctotal * ( 1+ KHCO3 / CH) ;
CCO3_inter = KHCO3 *CHCO3_inter /CH ;
Ca = Ca_total;
ph1 = CH + 2*Ca -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter + Kw/CH );
end
function y = SO2_OdeDriver (y0)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 KSO2 ...
KHSO3 Kw MWCaSO3 rhoCaSO3 P CCO2_in S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 CH ...
nloop Delt t CH_trial Tag pHInter
% take nloop steps starting at y0.
for k = 1:nloop
pH_trial = 8.0;
pH = fzero (@HionpH, pH_trial );
CH = 10^(3-pH);
% report as pH for ease of reference
pH =3 -log10 (CH) ;
% calculate the specieswise bulk concentrations by calling
%% C_Bulk
% Relationship between CSO2_bulkslurry, S_total,CH+,KSO2 and KHSO3-
CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between HSO3-, S_total,CH+,KSO2 and KHSO3-
CHSO3 = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between SO32-, S_total,CH+,KSO2 and KHSO3-
CSO3 = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between CCO2_bulkslurry, C_total,CH+,KCO2 and KHCO3-
CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);
% Relationship between HCO3-, C_total,CH+,KCO2 and KHCO3-
CHCO3 = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);
% Relationship between CO32-, C_total,CH+,KCO2 and KHCO3-
CCO3 = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
%%
% use these to find dydt expressions calling the odes.m file.
dydt = odes (t, y0);
% integrate with explicit Euler
y = y0 + dydt' * Delt;
% new results are
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);
% reassign these for the next time step.
y0 = y ;
t = t +Delt;
CH_trial = CH;
end
end
function dcdt = odes (t, c)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 KSO2 ...
KHSO3 Kw MWCaSO3 rhoCaSO3 CSO2_bulkslurry CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 CH CCO2_in S_total C_total Ca_total Trial ...
CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter CCO3_inter pHtrial pi pHInter pg
%% variable definitions for ease of reference
% c(1) = SO2 in the exit gas.
% c(2) = CO2 in the exit gas.
% c(3) = total sulfur in bulk liquid
% c(4) = total carbon in bulk liquid
% c(5) = total calcium in bulk liquid
% c(6) = solid CaCO3
% c(7) = solid CaSO3.
% SO2 in gas phase.
S_total = c(3) ;
C_total = c(4) ;
Ca_total = c(5) ;
% trial value
pH1 = 7;
pH = fzero (@HionpH, pH1 );
CH = 10^(3-pH) ;
%% C_Bulk
% Relationship between CSO2_bulkslurry, S_total,CH+,KSO2 and KHSO3-
CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between HSO3-, S_total,CH+,KSO2 and KHSO3-
CHSO3 = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between SO32-, S_total,CH+,KSO2 and KHSO3-
CSO3 = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);
% Relationship between CCO2_bulkslurry, C_total,CH+,KCO2 and KHCO3-
CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);
% Relationship between HCO3-, C_total,CH+,KCO2 and KHCO3-
CHCO3 = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);
% Relationship between CO32-, C_total,CH+,KCO2 and KHCO3-
CCO3 = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
%%
CH;
pH;
% SO2 in gas phase
% bulk gas at exit conditions
p_exit = c(1)*R*T ;
pg = p_exit;
Trial ;
Tnew = fsolve ( @NSub , Trial);
% Results are
pHInter = 3-log10(Tnew(3) );
pstar = Tnew(1) * HSO2;
NSO2 = kga *(pg-pstar);
% assign new trial values
Trial = Tnew;
dcdt(1) = (1/V_Headspace) * (F * CSO2_in - F * c(1) ) - NSO2 ;
% CO2 in gas phase.
E_CO2 = 1;
NCO2 = 0;
dcdt(2) = (1/V_Headspace) * (F * CCO2_in - F * c(2)) - NCO2 ;
% total sulfur balance
% CaSO3 cystallization sub-model
n = 3;
RCaSO3 = 0 * 1.62E-08 * exp(-5152/T) *( (c(5)*CSO3/ KSPCaSO3) -1)^n ;
% Limestone dissolution submodel
RCaCO3 = ktot*BETCaCO3*MWCaCO3*CCaCO3*CH *(1 - (Kad*CH)/(1+Kad*CH));
% RCaCO3 = NCO2 + RCaCO3;
dcdt(3) = NSO2 - RCaSO3;
dcdt(4) = NCO2+ RCaCO3;
% total calcium balance in liquud
dcdt(5) = RCaCO3- RCaSO3;
if ( c(6) < 1.0e-3 )
c(6) = 0.0; % cut off caco3 here.
dcdt(6) = 0.0;
Time_c = t ; % time when caco3 becomes zero.
Tag = 1;
else
dcdt(6) = - RCaCO3 ;
end
% CaSO3 balance
dcdt(7) = RCaSO3;
end
function discrep = NSub (Tr)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 KSO2 ...
KHSO3 Kw MWCaSO3 rhoCaSO3 S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter ...
CCO3_inter CCO2_bulkslurry CHCO3 CCO3 pg pHInter
Tr;
CSO2_inter = Tr(1) ;
pstar = CSO2_inter * HSO2 ;
S1Total = Tr(2) ; % CHSO3 + CSO3 at interface
CH = Tr(3) ; % h-ion concentration at interface
S_total = CHSO3 + CSO3 + CSO2_bulkslurry ; % in bulk for rate-Liq calculation.
CHSO3_inter = KSO2*CSO2_inter / CH ;
CSO3_inter = KHSO3 *CHSO3_inter /CH;
Ctotal = CHCO3 + CCO3 ;
CHCO3_inter = Ctotal * ( 1+ KHCO3 / CH) ;
CCO3_inter = KHCO3 *CHCO3_inter /CH ;
S_total_inter = CSO2_inter + CHSO3_inter + CSO3_inter;
df = S_total_inter- S_total ; % driving force due to reactions
kga ;
pstar;
pg;
Rate_gas = kga* (pg-pstar) ;
Rate_liq = kLa * df;
% discrepancies to be resolved.
discrep(1) = S1Total - (CHSO3_inter + CSO3_inter) ;
discrep(2) = Rate_gas - Rate_liq ;
% minus of H conce calculated
mCH_cal = 2*Ca_total -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter + Kw/CH ) ;
discrep(3) = CH + mCH_cal ;
end

  댓글 수: 11

표시 이전 댓글 수: 8
I recommend you add the 'Display', 'none' option combination to your fsolve() call.
After reftifying "\n", the error is gone, however, I only get one value per species.
I have used :
options = optimoptions('fsolve','Display','none');
Tnew = fsolve ( @NSub , Trial,options);
But it did not make a difference.

로그인 to comment.

태그

답변 수: 0


Translated by