{Coffee beans are roasted at a non-air flow oven situated at a roasting temperature Ta. Beans are introduced to preheated oven, and heat is transferred to the beans through radiation and natural convection. During roasting moisture is removed from bean surface. A single bean is considered, and the shape is assumed an equivalent sphere (of same volume). Transfer of heat and moisture is assumed to occur only in radial direction, along bean radius Rd. Size and thermal properties of the bean are given as constants and time-moisture dependent variables. The properties of roasting air are given in polinomial functions, obtained by plotting reference data. Time unit used here is minute. In this study, the surrounding air is assumed unaffected by the roasting, as such the bulk temperature and moisture content are constant. Thus, only the bean is considered, and a 1D sphere coordinate was chosen. The heat and mass fluxes are at the boundary, thus appear as POINT LOAD functions at the surface point. The boundary heat flux consists of radiative and (natural) convective heat flux, as well as the heat used for moisture evaporation. The evaporation starts when the temperature reaches 100 celsius degree, thus the latent heat of evaporation is given by the function SWAGE. The boundary mass flux is given by the moisture removal governed by a mass transfer coefficient, which is governed by the Sherwood number. Heat source term appears as "qreact", but is currently neglected to make the model simpler} TITLE 'Coffee roasting in non-air flow oven' COORDINATES sphere1 VARIABLES Temp {bean temperature, in degC} MC {bean moisture content, in kg H2O/kg solid} SELECT THREADS = 3 DEFINITIONS g = 9.81 {m/s2} {Bean properties} Rd = 0.00318 {bean equivalent radius, m} A = 4*pi*(Rd^2) {bean surface area, in m2} V = (4/3)*pi*(Rd^3) {bean volume, in m3} dens = 1169.8289 {bean solid density, kg/m3} Cp = 2410 {bean specific heat, J/kg K} Tso = 15 {initial surface temperature, degC} c = 3.3336 b = 4.278 keff = c*(1e-2)*dens*(MC/(1+MC))+b {J/m K min} {Air properties} Ta = 220 {uniform temperature of oven air and radiative surface, in degC} MCa = 0 {moisture content in air, in kg H2O/kg air} rhoa = (-5.10319e-14)*Ta^5+7.09789e-11*Ta^4 -4.2785e-8*Ta^3 + 1.5971e-5*Ta^2 - 0.00469815*Ta + 1.29187 {air density, kg/m3} Cpa = (-1.38794e-13)*(Ta^6) + (1.88153e-10)*(Ta^5) - (9.73596e-8)*(Ta^4) + (2.31691e-5)*(Ta^3) - (1.97287e-3)*(Ta^2) + (0.0754189*Ta ) +1005.88 {air specific heat, J/kg K} ka = -0.0000013567782*Ta^2 + 0.004516467*Ta + 1.418589 {air conductivity, J/m K min} Mua = -2.00886e-15*Ta^4 + 2.833104e-12*Ta^3 - 2.362914e-9*Ta^2 + 2.91744e-6*Ta + 0.001037736 {dynamic viscosity of air, kg/m min} Pr = Mua*Cpa/ka {Prandtl number} {Convective heat transfer} beta = 2/(Ta+Temp) Grheat = g*beta*(Ta-Temp)*(Rd^3)/((Mua/(60*rhoa))^2) {Grashof number for heat transfer} Nu = 2 + 0.43*((Grheat*Pr)^0.25) {Nusselt number} ha = Nu*ka/(V/A) {convective heat transfer in air, in J/m2 K min} Bi = (ha/keff)*(V/A) {Biot number} hconv = ha/(1+0.3*Bi) {convective heat transfer from air to the bean, in J/m2 K min} {Radiative heat transfer} sigma = 5.67037e-8 {Stefan-Boltzmann constant, in W/m2 K4} epsb = 0.98 {bean emissivity} hrad = 60*4*sigma*epsb*((0.625*(Ta+273.15)+0.375*(Tso+273.15))^3) {radiative heat transfer coefficient, in J/m2 K min} qreact = 0 Deff = -6.6114E-04*Ta + 7.416E-02 {in m2/min} Lvap = SWAGE(100-Temp, 2.79e6, 0, 500) {latent heat of evaporation, in J/kg} {moisture mass transfer coefficient} Dwa = 1e-3 Sc = Mua/(rhoa*Dwa) {Schmidt number} betam = 2/((dens*MC)+(rhoa*MCa)) Grmass = g*betam*((dens*MC)-(rhoa*MCa))*(Rd^3)/((Mua/(60*rhoa))^2) {Grashof number for mass transfer} Sh = 2 + 0.59*((Grmass*Sc)^0.25) {Sherwood number} kmass = Sh*Dwa*(V/A) {mass transfer coefficient for moisture, from bean surface to the air} INITIAL VALUES Temp = 15 MC = 0.12 EQUATIONS Temp : dt(Cp*dens*Temp) = DIV(keff*GRAD(Temp)) + dens*qreact MC : dt(MC) = dens*DIV(Deff*GRAD(Temp)) BOUNDARIES REGION 1 START(0) LINE TO (Rd) POINT LOAD(Temp) = (A/V)*(hrad+hconv)*(Ta-Temp) + Lvap*kmass*(MCa-MC) {heat flux = radiation + convection + evaporative latent heat} POINT LOAD(MC) = (A/V)*(kmass)*(MCa-MC) TIME 0 TO 10 by 0.0001 !MONITORS !PLOTS HISTORIES HISTORY(Temp) AT (0), (Rd) HISTORY(MC) AT (0), (Rd) END