clear all
clc
%% PROPERTIES
cps = 0.83; 
cpg = 1.1; 
cpa = 1.06; 
cpash=1; 
rhog=1.1; 
rhoa=1.1; 
hevap=2256.4; 
hreac=400; 
Mg=30; 
Ma=29;
TR = 480; 
Tok = 20; 
Toilc=60; 

%% INPUT
mbin = 1000; % input biomass flow kg/hod
wbinsug = 40:1:60; % biomass moisture % weight
for n=1:length(wbinsug)
    wbin=wbinsug(n);
%% DRYING
wbout = 5; % required output biomass moisture %
mwout = mbin*(wbin-wbout)/100; %  water flow output kg/hod
mbdry = mbin-mwout; % biomass flow output kg/hod

%% REACTOR - comparative model (function reakce7.m)
M0=[mbdry 0 0 0];
time=0:0.001:3;
[t,x]=ode23s(@reakce7,time,M0);
x1m=100.*x./mbdry;
figure (1)
plot(time,x1m)
hold on
legend('biomass','char','tar','gas')
grid on
title('Comparative model')
close 
tR=2.5;
tRindex=tR/(time(2)-time(1));

mbr=x(tRindex,1);
mcr=x(tRindex,2);
mlr=x(tRindex,3);
mgr=x(tRindex,4);

cpb=cpbmean((TR+Tok)/2);
Ereq=(mbdry/3600*cpb*(TR-Tok)+mbdry*wbout/100/3600*hevap+hreac*(mbdry*(1-wbout/100)-mbr)/3600); 

%% COMBUSTION reaction 1
LHVg=9.49*1000; 
cch4=0.08;
cco2=0.2;
cco=0.2;
ch2=0.38;
cn2=0.14;

Mch4=16.04;
Mco2=44.01;
Mco=28.01;
Mh2=2.02;
Mn2=28;

Mg=cch4*Mch4+cco2*Mco2+cco*Mco+ch2*Mh2+cn2*Mn2;
mg=101.325*1*Mg/8.314/273;
LHVg=LHVg/mg;

LHVb=15.6*1000; 

LHVc=27.5*1000;

EcR1=LHVc*mcr/3600; % Energy from char combustion kW
EgR1=LHVg*mgr/3600; %Energy from pyro gas combustion kW
EbR1=LHVb*mbr/3600; %Energy from biomass combustion kW

Egain=EcR1+EgR1+EbR1; %kW Energy from combustion
EbalanceR1 = EcR1+EgR1+EbR1-Ereq; % MW energy balance - positive good, negative bad
OilYieldR1 = mlr/(mlr+mcr+mgr)*100; %weight yield of oil 

%% ENERGY COMBUSTION
Toc=850; %°C 

Vg=mgr/Mg; 
lambda=1.5;
no2=lambda*(cco*0.5+ch2*0.5+cch4*2)*Vg; 
mair=no2/21*100*Ma+6*(mcr+mbr); 

wbash=4; 
wcash=7; 

mash=mbr*wbash/100+mcr*wcash/100; 
mflue=mair+mgr+mbr+mcr-mash;


Qa=mair/3600*(Toc-Tok)*cpa; % kW required to heat air
cpc=cpcmean((Toc+TR)/2);
Qc=mcr/3600*(Toc-TR)*cpc; % kW required to heat char
Qg=mgr/3600*(Toc-Toilc)*cpg; % kW required to heat gas
cpb=cpbmean((Toc+TR)/2);
Qb=mbr/3600*(Toc-TR)*cpb; % kW required to heat biomass

Qash=mash/3600*(Toc-Tok)*cpash; %kW to be taken out in ash
Qflue=mflue/3600*(Toc-Tok)*cpa; % kW to be taken out in flue gas

Eremaining=Egain-Qc-Qg-Qa-Qb-Qflue-Qash; % kW available power in flue gas
Efluegascooler=Qash+Qflue; %kW in flue gas output

%% SAND COOLER
ms=Eremaining/cps/(Toc-TR); 
Qscon=Eremaining-Ereq;
Tscooler=Ereq/ms/cps+TR;

%% DRYER
TDB0 = 20; %°C
RH0=70; %%
p=101.325; %kpa
[h0,x0]=psych(TDB0,RH0,p);

TDB1sug=50:1:80;
for m=1:length(TDB1sug)
    TDB1=TDB1sug(m);
RH1=0:0.1:100;

for i=1:length(RH1)
     [h1t(i),x1t(i),w1t(i),p1t(i)]=psych(TDB1,RH1(i),p);  
end
j=1;
while abs(x0-x1t(j))>0.001
    j=j+1;
end
h1=h1t(j);
x1=x1t(j);

TDB2=0:0.1:TDB1;
RH2=70;

for k=1:length(TDB2)
    [h2t(k),x2t(k),w2t(k),p2t(k)]=psych(TDB2(k),RH2,p);  
end
l=1;
while abs(h1-h2t(l))>1
    l=l+1;
end

h2=h2t(l);
x2=x2t(l);

mdryair=mwout/(x2-x0); %kg/h vzduchu do susarny
Qdryer=mdryair/3600*(h1-h0); %kW potrebych k suseni

BalanceDryer(m)=Qflue-Qdryer;
end
result(n)=max(BalanceDryer);
end
o=1;
while result(o)>0
    o=o+1;
end
wbincritical=wbinsug(o)