Dear moderator and members,

I am now working on a modelling study of coffee roasting process inside a non-airflow oven. The case is as follows:

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.

Here I attach the code file.

I tried to write the code based on my understanding of the program, but I found some errors which are:

1. When the time step is made bigger, there will be floating point overflow error.

2. I tried to make the time step smaller, but an error occurred, saying there was a fractional power of negative number. I found it coming from the ((Pr*Grheat)^0.25). The term (Pr*Grheat) is a multiplication of Prandtl number and Grashof number for heat transfer. Prandtl number should be a positive constant in my case (I checked it), since all the air properties are given as functions of air temp, and the air temp is constant. If such error occurred, it should be the Grashof number that went negative. I don't understand why, since the Grashof number is given by the function Gr = g*beta*(Ta-Temp)*(Rd^3)/((Mua/(60*rhoa))^2), (NOTE: g = gravity in m/s2 and Mua = dynamic viscosity in kg/m s). Thus, the possible negative value could come only from (Ta-Temp). Ta is the air temperature, and Temp is the bean temperature. The Gr is actually meant to determine the convective heat transfer, which is activated only at the bean surface (point in 1D coordinate), and thus Temp is the temperature at the boundary. If the Gr number became negative, the possibility is because the Temp become bigger than Ta.

How come such thing happen? Are there other reasons my codes lead to errors?

Can define or activate the variables such as Gr, hrad, etc. only at the bean surface, which is the point Rd (bean radius) in 1D sphere coordinate?

3. I also used similar non-dimensional number for the mass transfer, and there is another term (Sc*Grmass) which I am afraid will potentially lead to similar error.

This is the first time I used the program, and I found it quite powerful. I really appreciate if I can get some help in my study and learning process in using this program.

Thank you and regards,

Daniel