Modelling of a single drum dryer

22 views (last 30 days)
Mladen Petrovic
Mladen Petrovic on 8 Oct 2020
Answered: CARLOS ALBERTO on 1 Dec 2022
Good evening,
I curently working on modelling of a single drum dryer. My problem is that i am not able to write a proper code. I am trying to solve simulatenous heat and mass balance using ode45 and the results that i getting doesn't make any sense.My code is based on research paper which you can find in the attachement. Below you can find my code.
What i am doing wrong and how i am supouse to solve heat and mass transfer for this problem?
Many many thanks in advance for helping out a struggling beginner!
%Drum properties
D=1.5; %Drum diameter, [m]
L=3; %Drum length, [m]
dx=0.05; %Wall thickness of the E15/30, [m]
kd=50; %Drum thermal conductivity, [W/m.K]
V_rc=10; %Drum rotation velocity, [sec-1/rotation]
%Product feed
ds_f=0.25; %Dry solid feed, [%]
mc_i=1-ds_f; %Initial mositure content (wet basis), [%]
mc_I=mc_i/(1-mc_i); %Initial mosisture content, [kg/kg dry basis]
D0=2.4e-9; %Moisture diffusivity at 25 oC,[m2/s]
Ea=2280e3; %Activation energy, [J/kg]
T_f=90; %Feed temperature,[oC]
c_ps=1.28e3; %Heat capacity feed, [J/kg.K]; SANNE(COSUN)
d0=0.002; %Initial feed thickeness, [m]; Assumed
ro_s=995; %Starch density, [kg/m3]; A.S.Mudjamar, page 651
ro_f=ro_s; %Feed density,[kg/m3]
k_s=0.481; %Starch thermal conductivity, [W/m.K]
%Basic calculations
S_d=pi*D*L; %Drum surface, [m2]
F_m=S_d*0.75*d0*ro_f; %Mass of feed on 3/4 of the drum, [kg]
r_t=V_rc*3/4; %Drum residence time, [s]
R2=D+d0+dx; %Drum external diameter, [m]
%Saturated steam
p_ss=12; %Steam absolute pressure, [bara]
T_ss=180; %Steam temperature, [oC]
%Dry product
mc_f=0.08; %Final moisture content, [%]
mc_F=mc_f/(1-mc_f); %Final moisture content, [kg/kg dry basis]
%Physical properties
c_pw=4.18e3; %Specific heat water, [J/kg.oC]
c_pv=1.87e3; %Specific heat water vapor, [J/kg.oC]
c_pa=1.01e3; %Heat content dry air, [J/kg]
delta_hv0=2257e3; %Heat of evaporation of water, [J/kg]
%Water evaporation
T_e=100; %Temperature of evaporation. [oC]
w_v=3216; %water vapor,[kg/h]
s_wv=w_v/S_d; %Specific water evaporation, [kg/m2.h]
M=18; %Water molar mass,[g/mol]
%Antoine constans
A=23.2906;
B=3882.07;
C=229.928;
%Ambient conditions
RH=0.20; %Relative humidity, [%]
Ta=25; %Ambient temperature, [oC]
pa=101325; %Ambient pressure, [Pa]
k_air=26.24e-3; %Air thermal conductivity for 25oC, [W/m.K]
R=8.314; %Universal gas constant, [J/mol.K]
ps_w=exp(A-B/(C+Ta)); %Saturated vapour pressure, [Pa]
pv_p=RH*ps_w; %Partial vapor pressure, [Pa]
Ra=286.9; %Individual gas constant air, [J/kg.K]
Rw=461.5; %Individual gas constant water vapor, [J/kg.K]
ro_a=(pa-pv_p)/(Ra*(273+Ta))+pv_p/(Rw*(273+Ta)); %Density, [kg/m3]
mc_a=1000*0.622*pv_p/(pa-pv_p); %Moisture content of air
%Total heat
Q_wv=(s_wv*c_pw*(T_ss-T_e)+s_dp*c_ps*(T_ss-T_e)+s_wv*delta_hv0)*S_d/3600/1000; %Energy demand for evaporation, [kW]
E_r=(Q_wv*1000/delta_hv0); %Rate of evaporation, [kg/s]
W_rem=r_t*E_r; %Water removed, [kg]
%Air
u_a=40000; %Air volume flow, [m3/h]
um_a=u_a*ro_a/3600; %Air mass flow, [kg/s]
um_a1=u_a*ro_a; %Air mass flow, [kg/h]
v_a=u_a/3600/(L*D); %Air velocity, [m/s]
w_total=um_a*(1-1/(1+mc_a/1000))*3600; %Total water, [kg/h]
ni_a=1.849e-5; %Table value,[Pa.s] (https://www.engineersedge.com/)
Dw_air=0.242*10^(-4); %Moisture diffusivity in air,[m2/s]; Engineering toolbox
Re=ro_a*v_a*D/ni_a; %Reynolds number
Pr=0.7926; %Prandtl number at 25 oC and 1 bara
Nu=0.027*Re^(0.805)*Pr^(1/3); %Nusselt number, correlation based on Re number
h_air=Nu*k_air/D;%Convective heat transfer coefficient, [W/m2.K]
Gm=(Dw_air/R/(Ta+273)/0.001)*log(1/(1-0.0313)); %Molar flux, [mol/m2.s]
Sc=ni_a/(ro_a*Dw_air);
Kg=0.281*Re^(-0.4)*Gm/(pa*Sc^0.56); %Mass transfer coefficient, [mol/m2.s.Pa]
Q_air=(c_pa*Ta*(um_a1-w_total)+w_total*(delta_hv0+c_pv*Ta))/3600/1000; %Energy level of air,[W]
%Vapor out drum dryer
D_a=um_a*3600-w_total; %Dry air, [kg/h]
Tot_f=D_a+w_total+w_v; %Total flow, [kg/h]
Q_total=Q_air+Q_wv; %Total energy level, [kW]
H=Q_total*3600/D_a; %Enthalpy,[kJ/kg dry air]
m_Ctot=(w_total+w_v)/D_a*1000; %Total moisture content, [g H2O/kg dry air]
pv_po=pa*m_Ctot/1000/(0.622+m_Ctot/1000); %Partial vapor pressure out, [Pa]
T_out=41.4; %Outlet air temperature, [oC]
ps_wo=exp(A-B/(C+T_out));
RH_out=pv_po/ps_wo; %Relative humidity of outlet air, [%]
D_p=B/(A-log(pv_po))-C; %Dew poind, [oC]
%Air velocity variation:
u_a0=0; %Initial air volume flow, [m3/h]
u_af=60000; %Final air volume flow, [m3/h]
u_avar=(u_a0:100:u_af);
for i=1:length(u_avar)
%Air
um_a(i)=u_avar(i)*ro_a/3600; %Air mass flow, [kg/s]
v_a(i)=u_avar(i)/3600/(L*D); %Air velocity, [m/s]
Re(i)=ro_a*v_a(i)*D/ni_a; %Reynolds number
Nu(i)=0.027*Re(i)^(0.805)*Pr^(1/3);
h_air(i)=Nu(i)*k_air/D;%Convective heat transfer coefficient, [W/m2.K]
w_total(i)=um_a(i)*(1-1/(1+mc_a/1000))*3600; %Total water, [kg/h]
Q_air(i)=(c_pa*Ta*(um_a(i)*3600-w_total(i))+w_total(i)*(delta_hv0+c_pv*Ta))/3600/1000; %Energy level of air,[W]
%Vapor out drum dryer
D_a(i)=um_a(i)*3600-w_total(i); %Dry air, [kg/h]
Tot_f(i)=D_a(i)+w_total(i)+w_v; %Total flow, [kg/h]
Q_total(i)=Q_air(i)+Q_wv; %Total energy level, [kW]
H(i)=Q_total(i)*3600/D_a(i); %Enthalpy,[kJ/kg dry air]
m_Ctot(i)=(w_total(i)+w_v)/D_a(i)*1000; %Total moisture content, [g H2O/kg dry air]
pv_po(i)=pa*m_Ctot(i)/1000/(0.622+m_Ctot(i)/1000); %Partial vapor pressure out, [Pa]
T_out=41.4; %Outlet air temperature, [oC]
ps_wo=exp(A-B/(C+T_out));
RH_out(i)=(pv_po(i)/ps_wo)*100; %Relative humidity of outlet air, [%]
D_p(i)=B/(A-log(pv_po(i)))-C; %Dew poind, [oC]
end
figure
title('Convective heat coef, Total energy vs. Outlet air velocity')
yyaxis left % y-axis on the left indicates Convective heat coef
plot(v_a,h_air)
xlabel('Outlet air velocity [m/s]')
ylabel('Convective heat coef [W/m2.K]')
yyaxis right % y-axis on the right indicates Total energy
plot(v_a,Q_total)
xlabel('Outlet air velocity [m/s]')
ylabel('Total energy [kW]')
hold on
figure
title('Dew point, Relative humidity vs. Outlet air velocity')
yyaxis left % y-axis on the left indicates Dew point
plot(v_a,D_p)
xlabel('Outlet air velocity [m/s]')
ylabel('Dew point [oC]')
yyaxis right % y-axis on the right indicates Relative humiduty
plot(v_a,RH_out)
xlabel('Outlet air velocity [m/s]')
ylabel('Relative humiduty [%]')
hold on
T1=table(v_a',Re');
teta=[0, 270];
X=[mc_I mc_F];
[teta,X]=ode45(@Massbalance4,teta,X)
figure
plot(teta,X(:,1))
hold on
teta=[0, 270];
T=[100 150]
[teta,T]=ode45(@Heatbalance4,teta,T)
figure
plot(teta,T(:,1))
hold off
  2 Comments
CARLOS ALBERTO
CARLOS ALBERTO on 17 Nov 2022
Hi, do you have the files? table.m, Massbalance4.m" and "Heatbalance4.m"?

Sign in to comment.

Answers (1)

CARLOS ALBERTO
CARLOS ALBERTO on 1 Dec 2022
Hi!!! I need your help please!!!

Categories

Find more on Thermal Analysis 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!