A HOMOGENIZATION APPROACH FOR THE ROASTING OF AN ARRAY OF COFFEE BEANS∗

While the processes underlying the roasting of a single coffee bean have been the focus of a number of recent studies, the more industrially relevant problem of roasting an array of coffee beans has not been well studied from a modeling standpoint. Starting with a microscale model for the heat and mass transfer processes within a single bean during roasting, we apply homogenization theory to upscale this model to an effective macroscale model for the roasting of an array of coffee beans. We then numerically simulate this effective model for two caricatures of roasting configurations which are of great importance to industrial scale coffee bean roasting: namely, drum roasters (where the beans are placed in a rotating drum) and fluidized bed roasters (where hot air is blown through the beans). The derivation of the homogenization problem has been carried out in a three-dimensional rectangular geometry. Simulations are presented both for simplified one-dimensional arrays of three-dimensional beans (as these are easier to visualize), as well as cross sections of full three-dimensional arrays of beans (for the sake of verification). We also verify our simulation results against experimental data from the literature. Among the findings is that increasing the air-to-bean volume fraction ratio decreases the drying time for the array of beans in a linear manner. We also find that, in the case of a fluidized bed, an increase in the hot air inflow velocity will decrease the drying time in a nonlinear manner, with diminishing returns observed beyond some point for large enough air inflow velocities.

1. Introduction. The coffee industry has grown to be worth over $100 billion worldwide, with over 2.25 billion cups of coffee being consumed each day and a myriad of related socioeconomic implications [30,41,49]. Although the methods of production of coffee have remained similar for decades, manufacturers are constantly looking to improve and optimize the process, requiring ongoing fundamental research [31,34,43]. The coffee production process starts with removing the pulp of the coffee cherry, leaving the green coffee bean, which is then dried [12]. These beans are then roasted, usually for between 8 and 15 minutes. Coffee bean roasting is a dry heat treatment process, in which the flavor, color, and both physical and chemical structures of coffee beans are modified [9,19,65]. During the roasting process, heat, water, and various gases are transported, chemical reactions occur, and the geometry and mechanical properties of the beans also evolve over time [15,29,37,55,60,61].
There have been a variety of modeling studies on various aspects of coffee manufacturing, such as the brewing process [39,46,47,48]. Although understanding the coffee roasting process is fundamental to manufacturers being able to produce a better product, much of the research into understanding the roasting process has involved using experimentation or statistical analysis to obtain simple empirical models [3,59,63]. However, there have been recent efforts in the literature to formulate mathematical models of the coffee bean roasting process [23,24]. A system of partial differential equations was derived in [22] to describe the evolution of heat and moisture within a bean, although any explicit phase changes between water and water vapor were neglected, and moisture was treated in the bulk. The authors in [26,27] derived a mathematical model using conservation equations, considering three distinct phases within a bean and including the effects of multiphase flow and water evaporation. This work was extended in [25], with a more involved model being proposed to better describe the phase change between water and water vapor. The model of [25] also accounts for sucrose, reducing sugars, and other organic compounds, which are vital to the quality of the final product.
While the processes underlying the roasting of a single bean have been the focus of recent modeling efforts [22,25,26,27], the roasting of an array of beans has not been well studied from a modeling standpoint. The roasting of beans is a many-body process with heat transfer occurring during bean-bean collisions, collisions with the walls of the roaster, and via the air flow. There are two common types of roasting processes used in industrial scale coffee roasting: drum roasting and fluidized bed roasting [4]. In drum roasting [1,44], the beans are placed in a rotating drum. In fluidized bed roasting [10,33,45], hot air is blown through the beans. Figure 1 highlights the geometric differences between these two types of roasters. Although there may be motion of the beans, due to either the inlet flow or motion of the roaster, for our study we assume that the array of beans is fixed and that the volume fractions for both configurations are similar in size and static in time. In the present paper, we shall focus on the roasting of a static array of beans, under the assumption that the airflow is the dominant heat transfer mechanism via conduction between air and beans or convection due to the applied inlet flow. We shall use the behavior on the microscale of the roasting of a single coffee bean to derive effective equations on the macroscale of the roasting of an array of beans by way of homogenization theory. The remainder of this paper is organized as follows.
We first extend the model of [26] and derive equations for the heat and mass transfer during the roasting of a single bean, the air surrounding it, and the flow of air through a roaster in section 2. This will constitute the microscale model. In section 3, we use homogenization theory to derive macroscopic effective equations for air flow and heat and mass transfer in an array of roasting beans, obtaining the upscaled model for an array of beans. We discuss the numerical approach for solving the cell and macroscopic problems in section 4. We next use the upscaled model to study specific roasting configurations in section 5, where we also give a qualitative comparison of our results with data from experiments and carry out a sensitivity analysis on relevant parameters. We discuss our findings and possible avenues for future work in section 6.
2. Model description. We model the roaster as containing a uniform array of beans with air in between. Through conservation laws, we shall formulate a system of partial differential equations for the heat and mass transfer within the roaster, as well as consider the boundary conditions at the air-bean interface. We first list several caveats.
Although real coffee beans are not perfect spheres, and have a degree of anisotropy [28], we consider the idealized situation where all beans are spherical, of the same radius r b , and that the bean material is isotropic. Furthermore, we assume that all beans are static (they do not move) and remain in a periodic array configuration for all time. Therefore, the drum roaster motion is neglected, whereas simulations will preserve lengthscales and hot air inflow properties of the drum roaster. For both the drum roaster and the fluidized bed roaster, we assume no motion of the beans with the flow. Rather, the beans remain fixed, with the hot air flowing through the array of beans.
We neglect energy exchange between beans and the energy exchange on the walls of the drum or the fluidized bed roaster. We permit convection and conduction of heat through bean-air interactions, which are likely the dominant heat transport mechanisms on the timescales we consider. Convection appears to be the primary heat transfer mechanism, due to the fairly rapid timescale on which the beans are heated, dried, and roasted. Once the surrounding air is hot enough, conduction due to air-bean interactions is also present, but it is hard to distinguish the influence of conduction from convection, as the former relies on the latter first transporting heat from the external environment. Bean-bean convection is likely negligible compared to these other heat transfer mechanisms, and so is neglected, although one could certainly generalize the models we study to account for this additional heat transfer mechanism.
2.1. Multiphase heat and mass transfer model. We consider a roasting chamber which is a composite medium Ω ⊂ R 3 that consists of a periodic array of beans surrounded by air, where the bean centers are located a distance δ apart. Here, is the characteristic length of the roasting chamber, and δ 1 is a small dimensionless parameter. Given r b as the characteristic radius of each bean, we must have δ ≥ 2r b and hence r b to permit δ 1 (the characteristic lengthscale of a bean must be much smaller than the characteristic lengthscale of the roaster). We define the air phase of the medium by Ω a , the bean phase by Ω b , such that Ω = Ω a ∪Ω b , and the interface between each bean and the air surrounding it by ∂Ω b . We definẽ x ∈ R 3 to be the spatial variable of the system. A full list of dimensional parameters is provided in Table 1.
Following from the modeling performed in [25,26], we consider a representative small volume within the bean with three phases present: solid, liquid, and gas. The solid phase consists of cellulose and other organic molecules, we assume that the liquid phase consists of only water, and the gas phase consists of only water vapor. We neglect the production of CO 2 and the presence of other gases, since vapor production  [25] and there is large uncertainty in model parameters and kinetics associated with the production of other gas species [25]. As we will be incorporating an airflow model for the exterior of the roasting beans, we omit the boundary conditions presented in [25,26]. We define the constant porosity, φ, to be the volume fraction that the liquid and gas phases occupy, and the saturation,S b , to be the volume fraction of liquid water divided by the total volume of liquid water and water vapor. Thus, the volume fractions of the solid, liquid, and gas phases are given by 1 − φ, φS b , and φ(1 −S b ), respectively. Throughout this paper, we refer to quantities within the air by subscript a and quantities within the bean by subscript b. We denote the vapor pressure and temperature within the air and bean byP a ,P b andT a ,T b , respectively. We assume that the mass flux of water is negligible due to the low permeability and small effective water diffusivity within the bean. Hence, the loss of liquid water is a function of the evaporation rateĨ v only, and the conservation equation is given by The kinetics of this evaporation rate is motivated by the Langmuir evaporation equation [40]. The equilibrium vapor pressureP * v , also known as the sorption isotherm, is defined as the product of the water activity and the steam table of pure water [25]: Here, σ is the initial water saturation of the bean, the parameters A 1 , A 2 , A 3 are linked to the steam table of pure water, and C 1 and C 2 are associated with the choice of sorption isotherm employed (cf. [20,25,52]); see Table 1. From the assumption that the surface of the bean is dry and that the air contains no water, we prescribe the saturation to be zero on the surface between the bean and the air. We note that this is not a boundary condition, asS b is governed by the ordinary differential equation (1). We assume that the mass flux of water vapor in both the air and bean is due to the pressure gradient of the vapor [50]. We also assume that water vapor is advected by the flow within the air (which we assume is independent of the vapor pressure in the air) and that the evaporation of water only takes place within the bean. Hence, using the ideal gas law to define the densities ρ vi =P i m v /RT i for i = a, b, we derive the equations governing the pressures P a and P b , which come from conservation of water vapor: As suggested in [25], we take the permeability of gas within the bean, k b , to be of similar magnitude to the longitudinal permeability of gas within wood (as studied in [13]), since both wood and coffee beans principally consist of cellulose and exhibit similar timescales in the diffusion of their gas species [2]. We assume the binary diffusivity of water vapor in air per unit pressure, D va , to be constant.
We now consider conservation of energy within the bean and air phases, relating the rate of change of enthalpy to divergence of the advective and diffusive heat fluxes. As described in [26], we need not consider the enthalpy of the gas phase within the bean, and we can neglect the advection of heat within the beans. Within the air outside of the beans, heat is transferred by diffusion and is advected by the airflow. Hence, we obtain the following conservation of energy equations: where λ vĨv is the energy required to evaporate the water. Note that terms involving u are responsible for convection of heat, through the flow problem forũ.
At each air-bean interface, we assume continuity of temperature and heat flux, as well as continuity of vapor pressure and vapor mass flux. This is a different boundary condition than the one used in [27], as we are now modeling an array of beans rather than just a single bean suspended in air. Hence, the boundary conditions within the domain Ω are given by (5) where n is the unit normal vector pointing from the bean to the air. The air-bean interface conditions on the temperature are responsible for the transfer of heat from the air to the beans by conduction, as they allow the two distinct heat conduction problems (for air and beans, respectively) to interact. The relevant dimensional parameters in this model are again listed in Table 1.
2.2. Flow model. As well as the heat and mass transfer model, we need to describe the flow of air within the roaster. We assume that the air is an incompressible Newtonian fluid which rapidly tends to its steady state (i.e., quasi-static), and that the air cannot flow through the surface of the bean.
Assuming an air temperature at 200 • C, based on the size of r b and U given in Table 1, the Reynolds number is Re = O(10). This is around the upper bound for the validity of Stokes flow, but well below large Reynolds numbers resulting in turbulent flows. For faster velocities U, full Navier-Stokes would be required to accurately model the laminar flow, but we feel as though Stokes flow is sufficient to understand the problem we consider in a qualitative sense. As such, we employ Stokes flow in constructing our effective model. Hence, the governing air flow equations are given by the Stokes equations Here, µ a is the viscosity of the air, which we assume to be constant, andP T is the total pressure of the mixture of gases in the air and is assumed to be independent ofP a . Due to the incompressibility assumption, the temperature and vapor pressure within the air do not influence the air flow, with the primary role of the airflow being to deliver heat to the beans via convection. Of course, one could consider compressible flows governed by Navier-Stokes or more complicated turbulence models. However, as a first theoretical study on multibean coffee roasting, we aim to keep the fluid mechanics problem as simple yet physically relevant as possible.

Nondimensionalization.
We nondimensionalize the multiphase heat and mass transfer model (1)-(5) and the flow model (6a)-(6c) using the scalings in Table 2. We focus on the characteristic timescale where vapor diffusion via Darcy's law occurs. To do this, we nondimensionalize time byt = where p 0 P is the steam table pressure at roasting temperature (T ∞ ) and all other dimensional quantities are listed in Table 1. Furthermore, we nondimensionalize the air temperature T a to be on the unit interval usingT a = T 0 + (T ∞ − T 0 )T a , where T 0 is the roaster room temperature and T ∞ is the temperature of the hot air which enters the roaster through the inlet; a similar nondimensionalization is performed for the bean temperature, T b .
Using notation similar to that in [27], we give the heat and mass transfer model for the time evolution of the primary solution variables S b , P b , T b , P a , and T a in Table 2 Description of variables and their nondimensionalizations as motivated by [26]. The dimensional variables are represented with a tilde and the dimensionless variables without.

Variable
Description Nondimensionalization Vapor partial pressure in airPa = p 0 PPa Air flow velocityũ = U u P T Total gas pressure in airP T = µaU δ 2 P T nondimensional form: The nondimensional evaporation rate, I v , and sorption isotherm, P * v (T b , S b ), are given (as in [25]) by The boundary conditions are given in dimensionless form by All relevant dimensionless groupings are listed in Table 3. The equations (7a)-(7c) governing the saturation, temperature, and vapor pressure inside the bean have been taken directly from [27]. However, as shown in Table 2, we have nondimensionalized the spatial variable with respect to the characteristic length of the roaster, , rather than the radius of an idealized spherical bean as in [26]. We nondimensionalize the total pressure byP T = (µ a U/δ 2 )P T as in [16]. This pressure scaling is chosen to balance the pressure gradient over the macroscale (the roaster) with viscous forces over the microscale (each single bean). The dimensionless Table 3 Description and typical values of dimensionless parameters used in the model.

Dimensionless
Relationship to Typical value parameter dimensional parameters of parameter flow equations are thus given by We also assume constant homogeneous initial conditions in both the air and bean phases, 3. Homogenization. Although we have formulated governing equations for the heat and mass transfer within a single bean and the air surrounding it, we wish to understand how the roasting process occurs within an array of beans within a roasting chamber. Hence, we upscale our model, going from the microscale model of an individual bean surrounded by air to a macroscale model for the entire airbean medium within a roasting chamber. The formulation of a macroscopic model is achieved via homogenization, the process of deriving macroscopic equations for systems with fine microscopic structures. For general texts on homogenization theory, see [6,18,35,42,53]. Homogenization theory has been applied to a variety of multiscale heat and mass transfer systems in order to obtain effective macroscale models [5,7,21,54,56,58,64]. While there are a number of references available, we shall most closely follow the framework of the recent work in [8,17,16,62].
We assume that the beans are arranged in a cubic lattice structure, where the centers of each bean are located strictly δ apart. We divide the air-bean medium into periodic cells, each containing a single bean surrounded by air. We introduce y = x/δ as the microscale variable, treating x and y as independent. We define ω to be the unit cell, and the air and bean phases of the unit cell by ω a and ω b , respectively, such that ω = ω a ∪ ω b . The air-bean interface is defined by ∂ω b and the boundary of the cell by ∂ω a . The macroscale and microscale problem geometry are shown in Figure 2.
As we have introduced two spatial scales, the spatial derivatives vary according to (12) ∇ We denote the volume of the air and bean phases in each cell by θ a and θ b , with θ = θ a + θ b denoting the total volume of a cell. We define each cell to be the unit cube [0, 1] 3 . Since ∂ω denotes the boundary of the cell, we have θ = |ω| = ω dω = 1. We We then write the bean volume fraction as simply θ b = 1 − θ a and need only to prescribe the value of θ a henceforth. As we do not incorporate bean-bean interactions, for spherical beans we require that θ a > 1 − π/6 ≈ 0.47639, where θ a = 1 − π/6 corresponds to the case where the surface of the bean touches the boundary of the cell (and hence the neighboring beans). Therefore, we treat θ a as a control parameter, which falls within the range 1 − π/6 < θ a < 1.
3.1. Homogenized flow model. We consider the air flow first of all as it is assumed to be independent of the heat and mass transfer problem. The approach for the flow problem will be the same as in [17], and hence we omit details and give the cell problem and results.
We proceed with the asymptotic expansion of the flow velocity and pressure in powers of δ: At leading order, the pressure is found to be independent of the microscopic variable y, i.e., P T (x, t). At the next order, O(1), we obtain a solution taking the form whereP T is a function which is constant in y, and the matrix function K and vector function Π satisfy the cell problem Integrating (14a) over the unit cell air domain ω a yields the homogenized Darcy relation T , u * is the volumetric average fluid velocity defined by u * = ωa u (0) (x, y, t) dω, and K is a matrix containing information on the microscopic fluid flow, To close the homogenized system, we consider O(δ) terms to find the macroscopic incompressibility equation Equation (16) is known as Darcy's law [36], which describes fluid flow through porous media, stating that the velocity of the fluid is proportional to its pressure gradient. Under the assumption of spherical beans, the isotropy of the cell problem implies that the matrix K given in (17) may be written as [17] K = κI, where κ is a scalar quantity. We refer to the function κ as the effective permeability, and we have where N is the space dimension of the cell problem. For all simulations we choose N = 3, although we take N = 2 when plotting solutions to cell problems in section 4 as this is easier to display graphically.

3.2.
Homogenized heat and mass transfer model. We now return to the heat and mass transfer model, employing the same homogenization techniques to formulate effective governing equations for the entire air-bean medium. Full details are given in Appendix A.
We consider the asymptotic expansions in powers of δ: , meaning that at leading order, the temperature and pressure in the bean are both independent of the microscopic variable y and are homogeneous in the air and bean phases. Similarly, we find that S At O(δ), due to linearity of the problem, search for a solution taking the form , which involves a temperature cell problem for the components of the vector functions Γ 1 a and Γ 1 b , namely Γ 1 aj and Γ 1 bj , respectively. The temperature cell problem is given by for j = 1, 2, 3. Similarly, the pressure cell problem for Γ 2 aj and Γ 2 bj is given by for j = 1, 2, 3.
We use O(δ 2 ) equations to obtain the effective equations; again, see Appendix A for details. We wish to obtain a homogenized closed system of equations for the temperature, vapor pressure, and saturation of the air-bean medium which is valid as Using the new scalings outlined in Table 4, we write the homogenized conservation equations as (24) Table 4 Parameters used in the system (24)- (31) which differ between the two roasting configurations. Parameter values valid for the fluidized bed roaster and the drum roaster are provided. Other parameters which are the same in both roasting configurations are given in Table 3. Differences in most model parameters are due to the value of , the lengthscale of the roasting environment. For the drum roaster, this is = 0.2m, whereas for the fluidized bed, this is = 1.0m, and so this will modify some parameters between the two roasting configurations by a factor of five.

Dimensionless
Fluidized Drum parameter bed roaster roaster θa 0.9 0.875 1.06 0.879 3009 2485 Here, the effective macroscopic temperature diffusivity coefficient is given by while the effective macroscopic vapor pressure diffusivity coefficient is given by As was true for the effective permeability, under the assumption of symmetry due to spherical beans, the diffusivitiesD 1 andD 2 reduce to scalar multiples of the identity matrix I, sayD j = D j I, where D j is a scalar function, for j = 1, 2. Due to this symmetry, we have that where λ k a and λ k b are constants, and N = 1, 2, 3. We may then write the scalar quantities D j as To complete our homogenized system, we have the conservation of water equation, which we write as The initial conditions of the system are given by (31) S * = 1, T * = T in , P * = P in at t = 0.
The full list of parameters of the homogenized system (24), (25), (30), (31) (along with their reference values) can be found in Table 4.
4. Numerical approach. In this section, we discuss the approach taken for numerical simulations of the flow, temperature, and pressure cell problems, which are in turn used in the heat and mass transfer model simulations. We carry out all our simulations using the finite-element software package COMSOL Multiphysics version 5.3. We define the cells to be the unit cube [0, 1] 3 , containing a spherical bean at their centre. As we do not incorporate bean-bean interactions, we impose that θ a > 1−π/6, where θ a = 1 − π/6 corresponds to the case where the surface of the bean touches the boundary of the cell corresponding to δ = 2r b .
We solve the flow cell problem (15) and calculate the effective permeability, κ, as defined in (19), for different values of the air fraction θ a . The results are shown in Figure 3a for the N = 3 dimension case to represent the cell structure for a three-dimensional bean. We observe that κ → ∞ as θ a → 1, which is also seen in [17]. We solve the temperature cell problem (22) with the typical parameter value B 1 = 0.7 for different values of θ a . The results in Figure 3b show that D 1 → B 1 as θ a → 1, as expected from (29). Similarly, we solve the pressure cell problem (23) with B 2 = 12865. Again, as expected from (29), we see in Figure 3c that D 2 → B 2 as θ a → 1. So as the air fraction tends to 1, the effective thermal and vapor pressure diffusivities tend to the relative thermal and vapor pressure diffusivities of the air phase.
Once we have solved the cell problems to determine the effective parameters κ, D 1 , and D 2 , along with the homogenized flow model (16) and (18), we are able to carry out simulations of the homogenized heat and mass transfer model (24), (25), (30), (31). However, we must first define the geometry and boundary conditions for this system. We predominantly consider a one-dimensional roaster, assuming hot, dry air flowing in from the left-hand side of the roaster and out through the right-hand side. Higherdimensional geometries, such as a three-dimensional roaster, are briefly considered in section 5.2; we will show here that a one-dimensional geometry faithfully represents more realistic three-dimensional geometries. Hence, our boundary conditions are given by From (16) and (18), we deduce that the air flow is constant. Hence, we set u * = 1, as we have already nondimensionalized with respect to the inlet air velocity U (which in our case is the air velocity throughout the roaster).

Comparisons of roasting environments.
Having outlined the homogenization approach in section 3 and the numerical solution of the cell problems and upscaled macroscopic problems in section 4, we are now prepared to obtain results relevant to industrial coffee roasting configurations. As mentioned in section 1, there are two primary roasting configurations: the drum roaster and the fluidized bed roaster. We shall apply the results of previous sections to determine how an array of coffee beans will roast in each configuration. In the case of a drum roaster, there exist experimental results [22], and we shall compare our simulations to these.

Simulation of a drum roaster and model validation.
We first carry out a simulation which is representative of a drum roaster, to qualitatively validate our model with the experimental data from [22], in which samples of 3 g of Arabica green coffee beans were roasted at 200 • C in a drum roaster (EXPO 500/E, STA impianti, Bologna, Italy) for 14 minutes. The moisture concentration was determined at the 2, 4, 6, 8, 10, and 14 minute marks. The air speed was measured to be 0.02 m/s. As in  Tables 3-4. We plot the dimensional temperatureT , saturationS, and pressureP .  Figure 4b) and experimental data from [22].
[26], we can relate the water saturation, S * , to the moisture concentration, c, by To compare this to the experiential data, we take a spatial average of the moisture concentration, which we denote by c av = 1 |Ω| Ω cdΩ. From this, the initial saturation is σ = 0.0826, given that the initial moisture concentration in the experiment was 5400 mol/m 3 .
To represent the drum roaster (EXPO 500/E, STA impianti, Bologna, Italy), we take our one-dimensional macroscopic geometry to be of length 0.2 m. We fit our  Tables 3-4. We plot the dimensional temperatureT , saturationS, and pressureP . model to the data with only one fitting parameter, θ a , the air fraction, as this is an unknown. The parameters changed for use in these simulations are shown in Tables 3-4. We assume that our geometry of a static array of beans which have air being blown through them approximates a drum roaster, in which beans move through relatively slow-moving air.
The results of the simulation are shown in Figure 4, where the dimensional temperature, saturation, and vapor pressure are plotted against the spatial dimension and time. Initially, the air-bean medium starts off cool, with the beans containing water. As hot air is blown in from the left-hand side, a sharp heat and drying front propagates through the roaster from left to right, until the whole medium is dry and at the air inlet temperature, a result of the assumption that the beans within the roaster were fixed in a static periodic lattice. The vapor pressure reaches its peak on the right-hand side of the roaster, just before the medium becomes fully dry and heated to the air inlet temperature, at which the pressure falls rapidly to zero.
In Figure 5, we plot the simulation of the average moisture concentration c av from our model (24), (25), (30), (31), as derived from Figure 4b, against the experimental data from [22]. We observe a good qualitative fit, which suggests that the model does indeed describe the roasting process well, especially as only one fitting parameter was used. Note that the air fraction used to fit our model to the data is θ a = 0.875. We selected this value in consultation with Jacobs Douwe Egberts (JDE), and in order to show agreement with the results of [22], which were obtained using a drum roaster. As conversations with JDE suggested that θ a should be less than that for a fluidized bed, we considered θ a < 0.9. A value of θ a = 0.875 provided the best fit to [22], without having to artificially adjust other parameters.

Simulation of a fluidized bed roaster.
We now consider the parameter regime in which the model (24), (25), (30), (31) represents the heat and mass transfer within a fluidized bed roaster. We assume that our static array of beans approximates the dynamic movement of beans within a fluidized bed roaster. For this set of simulations, we consider a one-dimensional macroscopic geometry of length 1 m and an air inlet velocity of 0.1 m/s. Parameters are again as listed in Tables 3-4.
We see the results in Figure 6, observing behavior similar to that in Figure 4. We note that the increased air flow speed causes the heat and drying front to propagate more quickly through the roaster, even though the roaster domain is five times larger than the drum roaster domain used in Figure 4. There is also a notable increase in saturation before the beans are dried and a significantly greater vapor pressure.

Verification with 3D simulations.
In the previous subsections, we considered simplified one-dimensional geometries and claimed that these were representative of the full three-dimensional roasting chamber. To verify this, we now consider a three-dimensional macroscopic geometry of length 1 m in the x-direction and 0.2 m × 0.2 m cross-section which models a fluidized bed configuration. In addition to conditions already discussed for the one-dimensional configurations, the three-dimensional geometry employs no-flux boundary conditions along the cross-section edges. The simulation was carried out in slightly lower resolution than the one-dimensional simulations, and we focus on the first half of the simulation time to display our results.
In order to generate plots, we fix a line at the center of the domain (y c = 0.1, z c = 0.1) and a line near the external boundary of the cross section (y s = 0.01, z s = 0.1). We plot the temperature profiles for each in Figure 7. We find that there is very little difference between the temperature profiles near the center of the roaster and the boundary of its cross-section (the maximal difference, found within the drying zone, is half a degree Celsius, whereas outside the drying zone the difference is negligible). This implies that a one-dimensional macroscopic geometry faithfully represents a more realistic three-dimensional version, while having the additional advantage of significant computational efficiency. Similar agreement is found for vapor and pressure plots, and similar results are obtained for the drum roaster configuration, so we do not show those plots here.
We find that the drying front propagates through the roaster in a planar configuration, with the heat and mass transfer uniform in the coordinates orthogonal to this drying front. Of course, one should note that this will change if the walls of the roasting chamber are cooled, but as previously mentioned, with the exception of the hot inlet air, we neglect external thermal transport (say, from the atmosphere outside the roaster) in this study. In future work, it may be prudent to consider conduction boundary conditions, rather than no-flux boundary conditions, for the thermal field at the exterior boundaries orthogonal to the direction of heat convection. Furthermore, inclusion of more realistic curved geometry, rather than rectangular geometry, may influence the propagation of the drying front within the drum roaster. Still, we feel our one-dimensional results are qualitatively useful, particularly as effective models capturing the averaged roasting behavior of the beans.

Sensitivity analysis.
The homogenized heat and mass transfer model (24), (25), (30), (31) is dependent on many parameters, some of which coffee manufacturers are able to control to some extent-in particular the air inlet velocity, U, and the air fraction, θ a . Two of the key properties of the model that manufacturers are interested in are the drying time and the maximum vapor pressure. Hence, it would be instructive to see how varying U and θ a affect these two properties within the The absolute difference between these two temperature profiles is also shown, with the maximal difference found to be only half a degree Celsius. Similar small differences are likewise noted for the vapor pressure and water saturation solutions at each location. fluidized bed roaster. We define the drying time to be the time it takes for the saturation at the right boundary to reach 1% of its initial value.
In Figure 8, we plot the drying time and the maximum vapor pressure as the air volume fraction, θ a , varies. Apart from parameters which are dependent on either θ a or U, the parameters used are those in Tables 3-4 (with the exception that 1 = 10 −2 was used to reduce computation time), with a one-dimensional geometry of length = 1 m and boundary and initial conditions as previously stated. We observe that as θ a increases, the drying time decreases in a linear fashion, and that the maximum vapor pressure increases at a growing rate. This can be explained due to the fact that a greater air fraction increases the rate of heat transfer, which in turn causes the medium to dry more quickly, and a larger build-up in vapor pressure. At the same time, as θ a increases, there is less coffee per unit cell, and hence less mass to dry. It is sensible, then, that the beans roast slightly faster in the fluidized bed than in the drum roaster, since there is 25% more bean volume per unit cell in the drum roaster configuration.
We show the effects of varying the hot air inflow velocity, U, in Figure 9. We see that as U increases, the drying time decreases rapidly, and that the maximum vapor pressure increases at a steady rate. Again, these can be explained by an increased (a) Total drying time (b) Maximum vapor pressure Fig. 9. Sensitivity of total drying time and maximum vapor pressure to the air inlet velocity, U , for the fluidized bed roaster. We carry out a simulation using a one-dimensional macroscopic geometry of length = 1m, with air volume fraction, θa = 0.9, and other parameters as in Tables  3-4, except 1 = 10 −2 .
air flow driving the heat transfer within the medium. For larger velocities, there are diminishing returns, and this is perhaps due to the decrease in residence time, t ∼ θa U , since for large U, the residence time may no longer be rate-limiting relative to the heat and mass transfer processes occurring within the bean. Additionally, the assumption of a Stokes flow is most reasonable for U = O(10 −1 ). Near and above U = O(1), the Reynolds number, Re, becomes too large to justify the assumption of Stokes flow, with full Navier-Stokes required for accurate resolution of the fluid problem once Re = O(10 2 ). While U = O(10 −1 ) is close to operational inflow velocities in real roasters and is within the regime where the assumption of Stokes flow is reasonable, one should be careful when extrapolating our results for U = O(1) in the design of new roasters if the goal is to use faster inflow velocities to more rapidly roast coffee beans. Still, for U = O(10 −1 ), there are very clear improvements in drying time given a relatively small increase in the inflow velocity U.
It is natural to wonder if there is a correspondence between U and θ a , since for large U one might expect that the beans will move in some manner. In order to perform the homogenization procedure, we have assumed that the array of beans is static and hence cannot move with the flow. As such, the individual beans remain fixed in location, and the value of θ a is similarly fixed, regardless of the velocity. As such, U and θ a are independent parameters in our study. If the beans were allowed to move in a nonuniform manner (which is more realistic, but beyond what our homogenization approach permits, due to the need for spatial periodicity), then U may indeed have an influence over θ a and would further result in θ a becoming a function of time. This may be a topic to consider in future work, but as it would require direct simulation of a many-body problem for the beans, rather than homogenization, we do not explore this point further here.

Qualitative behaviors and dominant balances.
Having carried out a variety of simulations for both roasting configurations, we note that both temperature and water saturation vary rapidly near a localized drying front, which appears to correspond to the convection of enough hot air over the array of beans. On the other hand, the vapor pressure changes far more gradually, and globally rather than locally. To better appreciate why these behaviors are observed, we return to the effective equations (24)- (25).
As can be seen in Table 4, the dominant parameter groups include H 1 , H 3 , H 4 , H 5 = O(10 3 ). Additionally, observe that the effective diffusivity for the energy equation (24) is D 1 = O(1). This suggests that the effective thermal Péclet number, Pe eff , satisfies Pe eff = H4 D1 = O(10 3 ). Therefore, the dominant transport mechanism will be convection with the flow of hot air, with conduction due to diffusion highly localized within the narrow drying region. Away from the drying region, there is negligible change in S * and negligible influence from the diffusivity; hence the air temperature will at leading order be governed by a transport equation of the form (35) ∂T * ∂t + η∇ x · (u * T * ) = 0, where η = H4(1+A4S * ) H1+H3T S * . In the case of a fluidized bed, for small time and before the drying front arrives, S * ≈ 1, and therefore η = H4(1+A4) H1+H3T = 0.980. After the drying front passes, S * ≈ 0, and η = H4 H1 = 0.648. As such, there is more rapid convective transport of heat into the roaster ahead of the drying region. Behind the drying region, the temperature of the effective air-bean medium rapidly tends to the steady state temperature corresponding to the inlet boundary condition. Upon reaching this thermal equilibrium, the terms in the transport equation tend to zero, and heat remains uniform (until the roasting apparatus is turned off). The drying region in the effective model appears then as a narrow boundary layer between the cool and hot beans, behaving in a manner akin to the drying front previously studied for a single isotropic bean [27], with the primary difference being that this narrow layer acts on the effective air-bean medium. In the case of a drum roaster, we have η = 0.157 ahead of the drying region, and η = 0.104 behind it, so the relative qualitative discussion carries over from the fluidized bed, with the difference being that transport is slower for the drum roaster. That said, the drum roaster is smaller by a factor of five, and so these numbers are roughly in alignment with the fact that both roasters will roast their respective beans in a similar amount of time. Table 4 imply Pe eff = 4510 for the fluidized bed, whereas for the drum roaster we find Pe eff = 902. Therefore, while thermal convection is dominant in both configurations, it is even stronger in the fluidized bed. Indeed, when comparing the narrow drying regions in Figures 4a and 6a, we see that the transition region from cold to hot beans is even more narrow for the fluidized bed, due to the relatively larger effective thermal Péclet number.

Parameter values in
On the other hand, for the effective vapor pressure equation (25), we see that the effective diffusivity parameter satisfies D 2 = O(10 4 ), which is by far the most dominant term. The nearest term which balances this will depend on the roaster configurations. In the case of a fluidized bed, H 7 = 29.4, and this results in a dominant balance Therefore, diffusion of water vapor pressure is by far the dominant term, advection of water vapor is the second strongest term, and other terms are negligible. On the other hand, in the case of the drum roaster, H 7 = 5.87 and H 8 = 5.95 are of the same size, and hence the dominant balance in the effective vapor pressure equation (25) is with local saturation playing a larger role. For both cases, the diffusion term is still dominant, and this results in the more gradual change of water vapor in Figures 4c and 6c.
6. Discussion. We have formulated a heat and mass transfer model within a single bean and the air surrounding it, incorporating the single bean dynamics for an array of beans in a roaster. Using homogenization theory, we have derived a macroscopic model for the air-bean medium in order to model the effective properties of an array of coffee beans during the roasting process. We solved the system numerically, using parameter values which represent either a drum roaster or a fluidized bed roaster, and validated the model with experimental data for water concentration obtained in [22]. We find that temperature transport is primarily convective, except for within a narrow drying region where convection of hot air and air-bean conduction balance. Meanwhile, the build-up of water vapor pressure is dominated by diffusion, resulting in a more gradual build-up of pressure over time (until all beans are roasted, at which point the pressure is released).
Our results point to some key similarities and differences between the idealized drum and fluidized bed roasters we have considered and bring to light additional phenomena not seen in the single bean models. For both roasting configurations, there is a relatively sharp thermal gradient, with the temperature going rapidly from the initial (cool) temperature to the final (hot) temperature. Likewise, there is a relatively sharp drying front which propagates through the array of beans, analogous to the sharp drying front which was previously seen in the roasting of a single bean [27]. In both configurations, there is also a gradual vapor pressure build-up, which continues to grow until the beans reach terminal low water saturation. We note that this vapor pressure build-up is much less than what was found in the roasting problem for single beans with smaller gas permeability values [25,26]. Our configuration permits a larger gas permeability within the bean (which we take to be O(10 −14 ) m 2 ), as we model the air-bean interface more accurately than is possible when considering a single bean in isolation. The relatively small pressure build-up we observe is therefore consistent with maximal pressure predicted in Figure 4 of [25]. The sudden drop in pressure that we observe once the beans have dried is also consistent with Figure 3 of [25], and as discussed in [25], this is due to the inclusion of the sorption isotherm in our model.
A primary difference between the multibean roasting configurations is the time taken to dry all beans. Although the fluidized bed roaster is five times the size of the drum roaster, the array of coffee beans dries more quickly in the fluidized bed. This is seen when comparing Figures 4a and 6a, where the drum roaster has dried the beans within 800 s, yet the fluidized bed has dried the beans within 700 s, despite the much larger domain size. There are two sensible reasons for this. First, each unit cell in the drum roaster has 25% more bean volume than does a corresponding unit cell in the fluidized bed configuration; hence it takes slightly longer to dry this extra bean mass. Second, while convective heat transfer is dominant in both roasting configurations, it is largest within the fluidized bed, allowing for heat to more rapidly become available for the drying process.
Regarding industrial recommendations, we observe a decrease in the total drying time as the volume fraction of air increases, suggesting that configurations such as the fluidized bed will dry (and therefore roast) beans faster when there is sufficient hot air between beans. More generally, the lower the density of beans, the faster drying will occur. Of course, there is a trade-off here; one will want to roast as many beans as possible simultaneously, and due to the linear decrease of roasting time with air volume fraction seen in Figure 8a, finding a reasonable balance between the quantity of beans roasted and the time taken to dry all beans is seemingly feasible. There is a far more rapid decrease in roasting time with the inflow velocity of hot air, as seen in Figure 9a. This suggests that most of the decrease in roasting time happens before a velocity of around U = 0.2 m/s to U = 0.4 m/s, with only marginal decreases in roasting time after that. Therefore, expending effort on designing fluidized beds with much more rapid inflow velocities is likely to yield diminishing returns, beyond a point. That said, for larger values of U ∼ 1 m/s, we recall that the assumptions underpinning the choice of Stokes flow break down. Therefore, if higher velocities are to be considered in practice, they should be supported by simulations which utilize full Navier-Stokes flow valid for Reynolds numbers Re = O(10 2 ). Additionally, for faster flow, the assumption that beans remain stationary becomes less reasonable.
There are many additional features in the roasting process that have not yet been considered. For instance, in reality, beans are in constant motion, colliding with each other and the walls of the roaster chamber. Therefore, it may be more instructive to consider the spatial average values of the dependent variables in our model, as we did when comparing the simulations of our model to experimental data. We assumed that the temperature and vapor pressure within the beans and in the air were the same. Future models might relax this assumption, and thus we would be able to model a roaster where the temperature and vapor pressure within the air and beans influence each other, but evolve independently of one another. This would better describe the roasting environment, where the beans and air initially start off at different temperatures. If the flow is rapid enough so that the beans are suspended in air, then one could consider the problem of particles immersed in a flow. For finite domains, this may result in the formation of turbulent eddies, due to the geometric structure of the roaster. One might also consider more complicated fluid dynamics within the air-in particular the effect of turbulence on the heat and mass transfer within the roaster in the case where the fluid problem is in the high Reynolds number regime. Another option to obtain a more realistic model within a drum roaster, or other configurations where beans are in close proximity (or even touching), is to consider a bulk granular flow model. In addition to this, it would be beneficial to take into consideration the energy exchange on the walls of the roasting chamber.
We should finally note that similar modeling efforts to those applied to the roasting of a single coffee bean have been directed at understanding applications such as bread baking [66], carbon paste baking [57], wood drying [51], and popcorn popping [32]. While our motivating application has been to coffee bean roasting, the approach that we have taken may be adapted to the study of the roasting, baking, or drying of periodic arrays of material arising in these and other applications.
We consider the asymptotic expansions in powers of δ: We shall also use the following two expansions which will become relevant when expanding (38a) and (38b), respectively. First, (41) I and it suffices to state that Second, we make use of the expansion At O(1), we obtain the system , meaning that at leading order, the temperature and pressure in the bean are both independent of the microscopic variable y and are homogeneous in the air and bean phases. Similarly, from (45a), we can deduce that S b (x, t) due to the constant initial conditions (11).
At O(δ), we obtain Due to linearity of the problem, we may write the solutions as (47) T We can now formulate the temperature cell problem for the components of the vector functions Γ 1 a and Γ 1 b , namely Γ 1 aj and Γ 1 bj , respectively. The temperature cell problem then reads for j = 1, 2, 3. Similarly, the pressure cell problem for Γ 2 aj and Γ 2 bj is given by for j = 1, 2, 3.
From (47), we write The O(δ 2 ) terms in (38b)-(38e), (39a) become (52a) ∂ ∂t  = ∇ x · A 7 (∇ x T (0) a + ∇ y T (1) a ) + ∇ y · A 7 (∇ x T (1) a + ∇ y T (2) a ), y ∈ ω a , (52e) (∇ y T We now integrate (52a)-(52d) over the respective domains in which they are defined, the purpose being to use the divergence theorem to simplify the equations. We shall define dS as being the differential element of the bean surface ∂ω b . Any integrals evaluated on the boundary of the unit cell ∂ω a are zero, due to periodicity. First, integrating (52a) over ω b , and then using the divergence theorem and (50), we obtain (53) Similarly, integrating (52c) over ω a , using the divergence theorem and (50), we obtain The minus signs on the right-hand side arise because n is defined as the unit normal pointing outward from the bean to the air. The third and fourth terms on the lefthand side disappear as u (0) · n = u (1) · n = 0.
Integrating (52b) over ω b , using the divergence theorem and (50), gives (55) Using (43) in (46a), we have ∂S (1) b ∂t = 0, and from the initial condition (11) which gives S (1) b (x, y, 0) = 0, we deduce that S (1) b (x, y, t) = 0 for all t; hence the last term on the right-hand side of (55) is zero. Finally, integrating (52d) over ω a , using the divergence theorem, (50), and the volumetric average fluid velocity u * (as defined in section 3.1), we find We wish to obtain a homogenized closed system of equations for the temperature, vapor pressure, and saturation of the air-bean medium which is valid as δ → 0. Let us now denote T * = T (0) , P * = P (0) , and S * = S = P (0) . To obtain a system of homogenized equations, we must eliminate any higher order terms. We first formulate a homogenized conservation of energy equation. Taking the linear combination of A 7 multiplied by (55) added to A 3 B 1 (1 + A 4 S (0) b ) multiplied by (56), and using (52e), we eliminate higher order terms, and find (57) Here, the effective macroscopic temperature diffusivity coefficient is given by Similarly, to formulate a homogenized conservation of vapor pressure, we consider the linear combination of A 6 multiplied by (53) added to B 2 multiplied by (54), and using (46c), (46d), and (52f), we eliminate higher order terms, and find ∂ ∂t P * (A 6 θ b (1 − σS * ) + B 2 θ a ) 1 + T T * + B 2 A 5 ∇ x · u * P * 1 + T T * (59) where the effective macroscopic vapor pressure diffusivity coefficient is given by Using the new scalings outlined in Table 4, we simplify the homogenized equations, writing (57) as (61) [ and (59) as (62) ∂ ∂t To complete our homogenized system, we have the conservation of water equation, which we write as The full list of parameters of the homogenized system (61)-(63) (along with their reference values) can be found in Table 4. The initial conditions of the system are given by (64) S * = 1, T * = T in , P * = P in at t = 0.