Modelling of a single drum dryer

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

[Thanh D.B. Nguyen]

Hi,

I am very interested in the modeling of rotatry dryer.

I found that your above code needs three more additional matlab functions named "table.m", "Massbalance4.m" and "Heatbalance4.m".

Could you please send me these three functions then I can learn from your modeling experience.

Thank you very much in advance!

Thanh NDB,

[CARLOS ALBERTO]

Hi, do you have the files? table.m, Massbalance4.m" and "Heatbalance4.m"?

Answer 1

[CARLOS ALBERTO]

Hi!!! I need your help please!!!