Drying of Sisal Fiber: A Numerical Analysis by Finite-Volumes

: Vegetable ﬁbers have inspired studies in academia and industry, because of their good characteristics appropriated for many technological applications. Sisal ﬁbers (Agave sisalana variety), when extracted from the leaf, are wet and must be dried to reduce moisture content, minimizing deterioration and degradation for long time. The control of the drying process plays an important role to guarantee maximum quality of the ﬁbers related to mechanical strength and color. In this sense, this research aims to evaluate the drying of sisal ﬁbers in an oven with mechanical air circulation. For this purpose, a transient and 3D mathematical model has been developed to predict moisture removal and heating of a ﬁber porous bed, and drying experiments were carried out at different drying conditions. The advanced model considers bed porosity, ﬁber and bed moisture, simultaneous heat and mass transfer, and heat transport due to conduction, convection and evaporation. Simulated drying and heating curves and the hygroscopic equilibrium moisture content of the sisal ﬁbers are presented and compared with the experimental data, and good concordance was obtained. Results of moisture content and temperature distribution within the ﬁber porous bed are presented and discussed in details. It was observed that the moisture removal and temperature kinetics of the sisal ﬁbers were affected by the temperature and relative humidity of the drying air, being more accentuated at higher temperature and lower relative humidity, and the drying process occurred in a falling rate period.


Introduction
Vegetable fibers originate from plants and, basically, are composed of a solid skeleton and pores filled with fluid. The major components of these materials are cellulose, hemicellulose and lignin, with minor percentages of pectin, proteins, wax, inorganic salts and other water-soluble substances. Based on the cultivation region, soil type and climatic conditions, Energies 2021, 14, 2514 3 of 25 adequately controlling drying time and temperature. In this way, it is possible to obtain products with an acceptable quality level.
Diniz et al. [13], in their research related to the drying of sisal fibers, proposed lumped mathematical models to predict the transient behavior of the average moisture content and surface temperature, and also, the equilibrium moisture content of vegetable fiber as a function of the drying air temperature and water vapor concentration in the fiber porous bed. According to the authors, predicted results were in good fit with experimental data; for all drying conditions used correlation coefficients greater than 0.99 were obtained, and thus, was shown the potential of the proposed models to appropriately predict the drying process of sisal fiber, and other related materials.
Nordon and David [17] studied the coupled diffusion of heat and moisture in hygroscopic textile materials. The authors presented a 1D-numerical solution of the governing equations by finite differences, based on the "double scan" method for the nonlinear differential equations. The study was applied to wool, as an example of hygroscopic material and showed that the moisture transfers from the air to the wool and from the wool to the air are not symmetrical processes. According to the authors, the method does not predict the physical problem adequately.
Haghi [18] studied the simultaneous heat and moisture transfer in porous systems. For this purpose and based on the model reported by Nordon and David (1967) [17], the authors developed a one-dimensional mathematical model and its numerical solution using the finite difference method. The effects of operational parameters, such as temperature and moisture in the dryer, the initial moisture content of the porous material and the mass and heat transfer coefficients are examined using this model. The aim of this study was to describe adequately the heat and mass transfer during convective drying. According to the authors the proposed model can be used to predict transient variations in temperature distribution and moisture content in fabrics with reasonable accuracy.
Xiao et al. [19] and Xiao et al. [20] reported studies about mass transfer (imbibition process) in fibrous porous media using the fractal theory including surface roughness. The global model is based on the capillary theory (Hagen-Poiseuille model and Kozeny-Carman equation), assuming that the porous media is formed by a numbers of capillaries tubes. According to the authors, the proposed model predicts adequately the physical mechanisms of fluid transport inside fibrous porous media, and it contains no empirical constant, as commonly required for different models.
As already mentioned in previous comments, different distributed models have been reported in the literature applied to heat and/or mass transfer in fibrous porous media; however, restricted to one-dimensional approaches, no works used the finite volume technique for numerical solution of the governing equations. Then, in complement to these studies, this work aims to study the drying of sisal fibers in an oven using a 3D advanced mathematical modeling that considers the coupled phenomena of heat and mass transport inside the material, including fiber bed porosity and sorption heat. Here, the finite-volumes numerical method has been used to solve the governing equations. The idea is to assist engineers, industrials and academics in a rigorous understanding of the dominant mechanisms involved in the diffusion process occurring in the drying process. In this paper, the physical problem consists in predicting the moisture and heat transfer within a wet fibrous medium in the shape of a rectangular prism, as illustrated in the Figure 1. The fibrous medium is formed by a series of rigid and nonreactive fibers. In the drying process, hot air flows around the porous bed supplying heat to the wet fibers. The physical effect is to evaporate water inside the fiber, to heat the fibers, and to remove the water in the vapor phase. Water molecules in the vapor phase are free to move between Energies 2021, 14, 2514 4 of 25 the voids through the fiber bed, and to be absorbed or desorbed by the fibers. Volume variations provoked by water removal and heating of the fiber are not considered. the water in the vapor phase. Water molecules in the vapor phase are free to move between the voids through the fiber bed, and to be absorbed or desorbed by the fibers. Volume variations provoked by water removal and heating of the fiber are not considered.

The Governing Equations
Vapor Diffusion Equation To predict the water vapor concentration behavior inside the fiber bed along the drying process, the following equation was considered: where σ and β are constants from the equilibrium equation; ε is the porosity; ρ is the fiber density; C is the water vapor concentration (interfibers); T is the temperature, and t is the time. The parameter D = εD′ and D′ represents the vapor diffusion coefficient within the fibrous medium. In Equation (1), it is clear the effect of heat transfer in the mass transfer phenomenon. In the model, the following initial and boundary conditions were adopted: • Initial condition: • Boundary conditions: Heat Conduction Equation To predict the temperature of the fibers bed along the drying process, the following energy equation was used:

The Governing Equations Vapor Diffusion Equation
To predict the water vapor concentration behavior inside the fiber bed along the drying process, the following equation was considered: where σ and β are constants from the equilibrium equation; ε is the porosity; ρ s is the fiber density; C is the water vapor concentration (interfibers); T is the temperature, and t is the time. The parameter D = εD and D represents the vapor diffusion coefficient within the fibrous medium. In Equation (1), it is clear the effect of heat transfer in the mass transfer phenomenon. In the model, the following initial and boundary conditions were adopted: • Initial condition: • Boundary conditions: Heat Conduction Equation To predict the temperature of the fibers bed along the drying process, the following energy equation was used: where c p is the specific heat; K is the thermal conductivity; h s is the enthalpy and, T is the temperature. In Equation (6), it is considered that the energy of the wet air inside the fiber bed is negligible compared to the energy of the fibers. In Equation (6), the effect of mass transfer in the heat transfer phenomenon is clear.
• Initial condition: • Boundary conditions for heat transfer: Λ Equilibrium Equation According to Crank [21], the following equilibrium equation was considered: where M is the moisture absorbed per unit mass of fiber (kg vapor /kg dry fiber ), and α, σ and β are constants, that can be obtained by fitting to experimental data. From Equation (11), it is possible to obtain the following derivative: The average value of the potential of interest can be determined as follows: where Φ can be C, T or M (Equation (11)) and V is the volume of the fibrous medium.

Numerical Procedure
Equations (1) and (6) are highly nonlinear differential equations. Thus, exact solutions of these equation isn t possible. Thus, in this paper, these equations are solved using the finite-volume method [22][23][24]. In the discretization process, the following assumptions were adopted: (a) The solid is homogeneous and isotropic; (b) The only mechanism of water transport within the solid is diffusion; (c) The thermophysical parameters are constant along the drying process; (d) The diffusion coefficient is constant along the drying; (e) Convective mass and heat transfer coefficients are constant along the drying.
Due to the geometric shape (rectangular prism), it is possible to numerically solve the physical problem treated here using only one symmetrical part of the domain, as illustrated in Figure 2. In Figure 3 is shown the computational domain utilized in the discretization process of the governing equation. Based on the Figure 3, and integrating Equation (1) into volume and time, the following linear algebraic equation is obtained for the internal control volumes P: where,  For boundary volumes, the integration process of the governing equation is similar, as described for the internal points. However, must be considered the existing boundary conditions. For example, considering the symmetry condition, the mass fluxes in the west, back and south faces are null, i.e., C " = 0, C " = 0 and C " = 0, as illustrated in Figure 4. For the heat transfer equation after integration procedure, the following discretized linear algebraic equation is obtained: Based on the Figure 3, and integrating Equation (1) into volume and time, the following linear algebraic equation is obtained for the internal control volumes P: where, For boundary volumes, the integration process of the governing equation is similar, as described for the internal points. However, must be considered the existing boundary conditions. For example, considering the symmetry condition, the mass fluxes in the west, back and south faces are null, i.e., C w = 0, C b = 0 and C s = 0, as illustrated in Figure 4. For boundary volumes, the integration process of the governing equation is similar, as described for the internal points. However, must be considered the existing boundary conditions. For example, considering the symmetry condition, the mass fluxes in the west, back and south faces are null, i.e., C " = 0, C " = 0 and C " = 0, as illustrated in Figure 4. For the heat transfer equation after integration procedure, the following discretized linear algebraic equation is obtained:  For the heat transfer equation after integration procedure, the following discretized linear algebraic equation is obtained:

A T = A T + A T + A T + A T + A T + A T + B,
where, In the discretized form, Equation (6) applied to vapor concentration assumes the form: where, where i, j and k define the position nodal point, ∆V ijk = ∆x∆y∆z is the volume of the elemental volume, and npx-2, npy-2, and npz-2 are the number of control volumes in the x, y, and z directions, respectively. To solve iteratively (Gauss-Seidel method) the systems of algebraic equations originated from the Equations (14) and (23) as applied for all control-volume, a computational code using Mathematica ® software was developed. In the numerical treatment, the following convergence criterion was used at each nodal point and process time: where Φ can be C or T, and n represents the nth iteration at each time point. Simulations were performed using a numerical mesh of 20 × 20 × 20 nodal points and ∆t = 20s, which were obtained after a rigorous mesh and time step refining study.

Experimental Study
In this work, sisal fibers (Agave sisalana variety) were used in the drying experiments. The wet fibers were obtained in a farm located close to the town of Pocinhos, State of Paraiba, Brazil. The fibers were submitted to drying in oven with forced air circulation at different operating conditions. In Figure 5 is shown the sisal fibers bed used in the experiment. Table 1 summarizes information about the fiber and drying air.  In each drying experiments, the sample were withdraw from the oven periodically. In each step, the mass of was evaluated using a digital electronic device (0.001 g precision), and the surface temperature was measured using an infrared thermometer. These measurements were taken every 5 min (about 30 min), after 10, 15, 20, 25 and 30 min time intervals. After this period, the measurements were taken every 60 min until that constant mass was reached. Following that, the sample was dried for 24 h at the same drying temperature, in order to obtain the sample mass at the equilibrium condition, and then for a further 24 h at 105 °C to obtain the dried sample mass. Furthermore, at every measurement instant, temperature and relative humidity of the air outside of the oven were measured. This information was used to calculate air relative humidity inside the oven. In each drying experiments, air velocity inside the oven was measured using a hot wire anemometer. Figure 5 illustrates the point where surface temperatures were taking.

Transport Coefficients
In this research, the estimation of the mass diffusion coefficient and convective heat and mass transfer coefficients was performed using the least square error technique. Deviations between experimental and predicted values and variance were determined by [25]:   17 In each drying experiments, the sample were withdraw from the oven periodically. In each step, the mass of was evaluated using a digital electronic device (0.001 g precision), and the surface temperature was measured using an infrared thermometer. These measurements were taken every 5 min (about 30 min), after 10, 15, 20, 25 and 30 min time intervals. After this period, the measurements were taken every 60 min until that constant mass was reached. Following that, the sample was dried for 24 h at the same drying temperature, in order to obtain the sample mass at the equilibrium condition, and then for a further 24 h at 105 • C to obtain the dried sample mass. Furthermore, at every measurement instant, temperature and relative humidity of the air outside of the oven were measured. This information was used to calculate air relative humidity inside the oven. In each drying experiments, air velocity inside the oven was measured using a hot wire anemometer. Figure 5 illustrates the point where surface temperatures were taking.

Parameters Estimation 2.4.1. Transport Coefficients
In this research, the estimation of the mass diffusion coefficient and convective heat and mass transfer coefficients was performed using the least square error technique. Deviations between experimental and predicted values and variance were determined by [25]: where n is the number of experimental points,n is the number of fitted parameters, and Φ can be M or T. The initial guess of the convective heat transfer coefficient was determined by common correlations for the average Nusselt, Reynolds and Prandtl numbers [26], applied to air flowing over a plane plate by: where, Nu j = 0.664 × Re j , k is the air thermal conductivity and j = 1,2 and 3 ( Figure 2), valid in the intervals 5 × 10 5 < Re ≤ 1 × 10 8 .
The initial guess of the convective mass transfer coefficient applied to the air was determined by common correlations for average Sherwood numbers and Schmidt, by [26]: where, Sh j = 0.664 × Re j
With the values of T eq and M eq for each experimental condition, a linear fit of Equation (39) to the experimental data was performed by the quasi-Newton method using the Statistica ® Software (convergence criterion 9.9 × 10 −5 ). After this statistical procedure the values of the parameters a 1 and a 2 were obtained.
On the other hand, considering the drying air as an ideal gas, the air density inside the fibrous bed was calculated by: where P is the atmospheric pressure, ρ air is the air density, R is the particular gas constant (atmospheric air), and T abs is the absolute temperature of the drying-air. Then, from Equation (40), the equilibrium water vapor concentration between the fibers, can be determined by: where AH is the air absolute humidity on the voids of the porous bed.

Estimation of the Fiber and Sample Parameters
The apparent density of the fiber bed was calculated by: where V = 2R 1 × 2R 2 × 2R 3 is the volume of the fiber bed.
To obtain the bed porosity, the following equation was used: where V fiber = N × L × πd 2 4 is the total volume of fibers, L and d, are the length and diameter of the sisal fiber, respectively, and N represents the number of fibers ( Figure 6). where V = N × L × is the total volume of fibers, L and d, are the length and diameter of the sisal fiber, respectively, and N represents the number of fibers ( Figure 6). The parameters thermal conductivity and specific heat of the fiber bed were determined as follows: and The parameters thermal conductivity and specific heat of the fiber bed were determined as follows: and Table 2 summarizes the geometrical data of the fiber and samples before drying [11,[13][14][15][16].

Fiber Morphology
The scanning electron micrographs (SEM) of untreated sisal fibers in moist condition are illustrated in Figure 7. Upon analyzing this figure, it can be verified that the arrangement of this particular fiber is similar to other natural fibers. It can be seen a spongier aspect, voids, rougher surface and a thin and compacted cellular arrangement. Furthermore, it can be observed that parenchyma cells are widely distributed along the fiber. The fibers presented diameters almost constant along the length and, minor differences between the fibers mean diameter were observed.

Analysis of Equilibrium Equation
Equation (46) predict the linear dependence of the moisture content with the temperature and water vapor concentration, fundamental in predicting the coupling between the Equations (1) and (6). A correlation coefficient R = 0.99 was obtained after the performed linear regression. Most details about the statistical procedure can found in the literature [11,[13][14][15][16]. arrangement of this particular fiber is similar to other natural fibers. It can be seen a spongier aspect, voids, rougher surface and a thin and compacted cellular arrangement. Furthermore, it can be observed that parenchyma cells are widely distributed along the fiber. The fibers presented diameters almost constant along the length and, minor differences between the fibers mean diameter were observed.

Analysis of Equilibrium Equation
Equation (46) predict the linear dependence of the moisture content with the temperature and water vapor concentration, fundamental in predicting the coupling between the Equations (1) and (6). A correlation coefficient R = 0.99 was obtained after the performed linear regression. Most details about the statistical procedure can found in the literature [11,[13][14][15][16]. Figures 8 and 9 illustrate the predicted and experimental average moisture content (dry basis) of the fiber sample as a function of time for the drying air temperature of 50 °C and 80 °C. Upon analyzing these figures, it can be seen that there is a good concordance between the simulated and experimental of the average moisture content along the process for every condition. Besides, it was verified that the average moisture content decreases with time until to reach the hygroscopic equilibrium condition. Furthermore, the drying rate increased and equilibrium moisture content decreased when higher drying air temperature and lower air relative humidity were used in the experiment. Therefore, this physical situation resulted in a short drying time.

Drying and Heating Kinetics
In Figures 10 and 11 are shown the comparison between the predicted and experimental transient values of the temperature at the sample surface at 50 °C and 80 °C. Upon analyzing these figures, it can be noted that there is good concordance between surface temperature values in the first 4000 s of process followed by minor discrepancies for long time. For example, in the drying at 50 °C and 27,400 s elapsed time difference, between the predicted and experimental surface temperature data there is 2.8 °C  Figures 8 and 9 illustrate the predicted and experimental average moisture content (dry basis) of the fiber sample as a function of time for the drying air temperature of 50 • C and 80 • C. Upon analyzing these figures, it can be seen that there is a good concordance between the simulated and experimental of the average moisture content along the process for every condition. Besides, it was verified that the average moisture content decreases with time until to reach the hygroscopic equilibrium condition. Furthermore, the drying rate increased and equilibrium moisture content decreased when higher drying air temperature and lower air relative humidity were used in the experiment. Therefore, this physical situation resulted in a short drying time.     In Figures 10 and 11 are shown the comparison between the predicted and experimental transient values of the temperature at the sample surface at 50 • C and 80 • C. Upon analyzing these figures, it can be noted that there is good concordance between surface temperature values in the first 4000 s of process followed by minor discrepancies for long time. For example, in the drying at 50 • C and 27,400 s elapsed time difference, between the predicted and experimental surface temperature data there is 2.8 • C difference, while at 80 • C and 17,800 s elapsed time, the difference in these values is 4.5 • C. These discrepancies may be attributed to minor errors in the temperature measurement procedure in each the experiments or even by the accuracy in the device itself. Besides the differences verified for long time may also be associated with the consideration of constant convective heat coefficient and null dimensional variations along the process as established in the mathematical modeling.  Upon analyzing the Figures 7 and 8, it can be seen that the temperature at the surface of the fibrous bed changed along the drying, reaching the thermal equilibrium condition faster at the highest air temperature. The increases in the surface temperature of the fiber with progress of drying indicates the moisture removal occurred on the falling drying rate period, i.e., the water flux inside the fiber to its surface is less than the water flux removed by the heated air at the surface. Furthermore, it was observed that the fiber temperature increased as the drying rate decreased until the equilibrium condition. These phenomena were cited in the literature [10,11,[13][14][15][16].

Water Vapor Concentration and Temperature Distributions
In Figures 12 and 13 is illustrated the water vapor concentration field inside the fiber bed, at the plans x = 0.025 m (R1/2) and y = 0.0125 m (R2/2), at the moments 200 s, 700 s, 5000 s and 8000 s of the drying process, respectively.
By analyzing the Figures 12 and 13, it is possible to see that water vapor concentration presented the highest values in the center of the fiber bed at any time. Upon analyzing the Figures 7 and 8, it can be seen that the temperature at the surface of the fibrous bed changed along the drying, reaching the thermal equilibrium condition faster at the highest air temperature. The increases in the surface temperature of the fiber with progress of drying indicates the moisture removal occurred on the falling drying rate period, i.e., the water flux inside the fiber to its surface is less than the water flux removed by the heated air at the surface. Furthermore, it was observed that the fiber temperature increased as the drying rate decreased until the equilibrium condition. These phenomena were cited in the literature [10,11,[13][14][15][16].

Water Vapor Concentration and Temperature Distributions
In Figures 12 and 13 is illustrated the water vapor concentration field inside the fiber bed, at the plans x = 0.025 m (R 1 /2) and y = 0.0125 m (R 2 /2), at the moments 200 s, 700 s, 5000 s and 8000 s of the drying process, respectively.
By analyzing the Figures 12 and 13, it is possible to see that water vapor concentration presented the highest values in the center of the fiber bed at any time. Thus, moisture flows from the center to surface. Besides, vapor concentration decreased at any position and time, tending towards its equilibrium condition for long drying time. The regions close to the edge dry faster than the others regions inside the porous bed, especially in the vertex region. Figures 14 and 15 show the distribution of temperature inside the fiber bed, analyzed in the plans x = 0.025 m (R 1 /2) and y = 0.0125 m (R 2 /2), at the moments 200 s, 700 s, 5000 s and 8000 s of the drying process, respectively. at any position and time, tending towards its equilibrium condition for long drying time. The regions close to the edge dry faster than the others regions inside the porous bed, especially in the vertex region. Figures 14 and 15 show the distribution of temperature inside the fiber bed, analyzed in the plans x = 0.025 m (R1/2) and y = 0.0125 m (R2/2), at the moments 200 s, 700 s, 5000 s and 8000 s of the drying process, respectively.      Upon analyzing these Fs, it can be seen that the temperature has the lowest results in the center of the fibrous medium along the process. The temperature also increased in any position, tending to the thermal equilibrium condition for long drying time. This behavior proved that heat flux occurs from the surface to the center of the porous sample. Furthermore, it is observed that the maximum temperature occurred close to the edge, especially in the vertex of the sample.
From already mentioned results, the regions close to the vertex dry and heat faster. Therefore, these regions are most susceptible to higher moisture and temperature gradients, which originate hydric and thermal stresses. Based on the intensity of these stresses, the material can achieve unfavorable deformation and rupture levels, and poor quality for a specific application. Under the industrial aspect these results are important to assist engineers in making decision about how the material must be dried and the fiber distributed in the bed, in order to minimize these effects.

Estimation of Transport Coefficients (D, hm and hc)
During drying process, at the surface of the fiber sample different phenomena occur simultaneously: heat and mass convection due to hot air flowing over the sample, and heat and mass transfer due to the water evaporation. Analyzing the effects of convective heat and mass transfer at the surface of the porous sample, boundary conditions of third Upon analyzing these Fs, it can be seen that the temperature has the lowest results in the center of the fibrous medium along the process. The temperature also increased in any position, tending to the thermal equilibrium condition for long drying time. This behavior proved that heat flux occurs from the surface to the center of the porous sample. Furthermore, it is observed that the maximum temperature occurred close to the edge, especially in the vertex of the sample.
From already mentioned results, the regions close to the vertex dry and heat faster. Therefore, these regions are most susceptible to higher moisture and temperature gradients, which originate hydric and thermal stresses. Based on the intensity of these stresses, the material can achieve unfavorable deformation and rupture levels, and poor quality for a specific application. Under the industrial aspect these results are important to assist engineers in making decision about how the material must be dried and the fiber distributed in the bed, in order to minimize these effects.

Estimation of Transport Coefficients (D, h m and h c )
During drying process, at the surface of the fiber sample different phenomena occur simultaneously: heat and mass convection due to hot air flowing over the sample, and heat and mass transfer due to the water evaporation. Analyzing the effects of convective heat and mass transfer at the surface of the porous sample, boundary conditions of third kind were considered in this research. Thus, two phenomena occur, simultaneously: moisture flux at the surface that is proportional to the difference between the moisture content at the surface of the fiber sample and the equilibrium moisture content on the air temperature, and heat flux at the surface that is proportional to the difference between the temperature at the surface of the fiber sample and the drying air temperature. For these situations, the convective heat transfer coefficient (h c ) and convective mass transfer coefficient (h m ) were obtained using the least square error technique. In this technique, the predicted and experimental data of fiber average moisture content and surface temperature during drying at different conditions are compared. The best values of these parameters correspond to the lowest value of the variance statistical parameter.
As already described earlier, from these initial values of the convective heat and mass transfer coefficients, an iterative mathematical procedure was performed to determine the optimized values of these coefficients, including the mass diffusion coefficient. Tables 6  and 7 summarize the values of the parameters estimated, relative error and variance for all experimental tests obtained after the statistical procedure.  After analysis of the Tables 5 and 6, small errors and variances can be seen, proving that the all mathematical procedures used in estimating the transport coefficients are appropriate. Further, we verified that both convective heat and mass transfer coefficient have small values typical of natural convections, due to small values of the air velocity inside the oven, and tend to increase when air drying temperature is increased. Some oscillatory behavior was verified in the convective mass transfer coefficient that can be attributed, for example, to the fact that the effects of dimensional variations due to heating (dilation) and moisture removal (shrinkage) were not considered, and the assumptions of constant fiber bed porosity and equilibrium equation has linear behavior during the drying process.
Particularly, the mass diffusion coefficient (D) was obtained by multiplying the vapor diffusion coefficient within the fibrous bed by the bed porosity, as already mentioned before, i.e. D = εD'. Upon analyzing of the Table 6, we verified that major mass diffusion coefficients are found when higher air temperature was used. Furthermore, since the porosity of the medium presented minor variation in the experiments (0.915 < ε < 0.922), the drying of the fibrous bed presented similar behavior than the drying of individual fibers, which can be confirmed from the value of the parameter D which is less than the parameter D AB (Table 3) for all experimental tests, at least an order of magnitude. Therefore, we state that the water vapor flux within the fibrous bed in less than in the air outside the porous sample [27,28].
Finally, it is well known that the drying process requires high energy consumption and is responsible for high pollutant emissions in the atmosphere due to the use of different fuels for drying-air heating. Thus, it is important to control drying systems in laboratory or industrial scales, minimizing both energy consumption and environmental impact.
In this research, we develop macroscopic and advanced heat and mass transfer equations applicable for capillary-porous bodies with parallelepiped configuration. Both the model and the finite-volume method used herein demonstrated great potential, being accurate and efficient to be applied in many practical physical problems such as: a.
Diffusion processes of wetting, drying, heating and cooling, whether coupled or separate; b.
Coupled diffusion processes of liquid, vapor and heat in different porous bodies; c.
It is possible to study different diffusion problems assuming variable or constant thermophysical properties, and equilibrium or convective boundary conditions; d.
Under numerical aspects it is possible to use uniform or nonuniform grids, and to consider changes in the dimensions (dilation and shrinkage) of the bodies in the transient simulations. e.
The finite volume method is unconditionally stable even when applied to solving nonlinear partial differential equations; f.
No restriction about the nature of the porous body bed (fruits, grains, vegetables, textiles, etc.) is required when it is considered as continuum media; g.
Good estimation of process time to dry the product, contributing to reduction in energy consumption and increase in productivity of dried porous bodies; h.
Rigorous checkup and understanding of the effect of process variables on product quality during drying.
Despite the advantages presented above, some limitations of the modeling can be cited: (a) Chemical transformation, thermal and structure inside the porous body that, in general, are responsible to provoke no uniform distribution of void are not considered; (b) Inadequacy to be applied in drying processes where high air-drying temperature and effects of pressure are important; (c) Mass diffusion is the only moisture migration mechanism, and gravity, capillarity, and filtration effects, for example, are not considered; (d) It is necessary to use estimated parameters after fitting nonlinear regression or other techniques for this purpose.

Conclusions
From the predicted results, it can be concluded that: (a) The proposed mathematical model proved to be useful for describing the drying process of sisal fibers, considering the good agreement between the predicted and experimental data of average moisture content and surface temperature of the sisal fibers obtained in all drying conditions; (b) Despite the highly nonlinear character of the governing equations that represent the mathematical model, the finite-volume method proved to be useful to solve the cited equations and to assist in predicting the phenomenon of heat and mass transfer within the porous sample; (c) Drying of sisal fibers occurred at a falling drying rate period; (d) The largest water vapor concentration and temperature gradients are located in the regions near the vertex of the fiber porous bed. Then, these regions are more suitable to deformation and hydric and thermal effects, which are responsible for reducing product quality after drying process.
Finally, in continuation of this research, the authors strongly recommend studying drying problems occurring with dimensional variations of the bed and intermittent drying that is a technique useful in minimizing energy consumption.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available upon request from the authors.

Acknowledgments:
The authors are grateful to the Federal University of Campina Grande (Brazil) for the research infrastructure and the references cited in the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Average Prandtl number --R 1 Half-length of the bed in the x-direction m R 2 Half-length of the bed in the y-direction m R 3 Half-length of the bed in the z- Distance between nodal points in the x-direction m δy n Distance between nodal points in the y-direction m δy s Distance between nodal points in the y-direction m δz b

Abbreviations
Distance between nodal points in the z-direction m δz f Distance between nodal points in the z-direction m δx e Distance between nodal points in the x-direction m ∆t Time step s ∆V ijk Volume of the elemental volume, m 3 ∆x Length of the control-volume in the x-direction m ∆y Length of the control-volume in the y-direction m ∆z Length of the control-volume in the z-direction m ε Bed porosity µ Air viscosity N·s/m 2 ρ s Dry solid density kg/m 3 σ Constant of the equilibrium equation m 3 /kg dry fiber Φ Potential of interest -Φ Average value of the potential of interest --