How to use fprintf
3 views (last 30 days)
Show older comments
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 Comments
Answers (0)
See Also
Categories
Find more on Testing Frameworks in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!