A 3D transient CFD model to predict heat and moisture transfer in on-farm stored grain silo through parallel computing using compiler directives: Impact of discretization methods on solution efficacy

Abstract Computational fluid dynamics model parameters highly depend on the storage structure, grain conditions, and airflow properties. Solution methods obtained for the lab-scale bins cannot be extrapolated for a relatively large on-farm silo. To validate the above hypothesis, a model was developed and validated with stored barley aeration in a 1000t silo. A mathematical function was used to initialize the discrete initial conditions in the silo followed by using DEFINE_macros to execute parallel computing. Time-step analysis was conducted followed by optimization of the solution methods in terms of the accuracy and time frame of the model’s output. Results showed that with a time-step of 4 s, temperature and moisture error were 1.1–1.7 °C and 0.3% wb while mean relative deviation were 3.1–4.4% and 2.2%, respectively. COUPLED algorithm resulted in a similar accuracy as of SIMPLE scheme except within spatial and transient formulations. However, the former algorithm was observed to take almost 190% more time than the latter scheme, limiting the simulation efficiency in an on-farm silo. Green Gauss Node-based gradient technique was found to be the appropriate for discretizing large silos. Model showed 45 kJ kg−1 energy emission that decreased the cooling potential of air in silo. As field evaluation of aeration strategies are time-consuming, this model can be used to obtain results that could shape the stored grain management practice.


Introduction
On-farm grain aeration is a non-chemical technique indispensable in stored grain handling and management practices. It is accomplished by blowing ambient air through the grain using a fan attached to the storage silo plenum. This alters the grain ecosystem through appropriate heat and moisture transfer mechanisms which minimize the risk of insect infestation or favorable conditions for mold growth by cooling or drying the grain, respectively and ensure safe storability.
Stored grain researchers have long sought effective methods to model the stored grain ecosystem to accomplish difficult tasks such as exploring the thermo-physical phenomenon during aeration throughout the stored bulk. Recently, Panigrahi et al. [1] have outlined the benefits of different numerical methods to develop models (since the 1970s) that have significantly narrowed this issue, proliferating in-depth research in developing grain management strategies. However, when it comes to accurate prediction of temperature and moisture distribution (as a function of velocity vector) in the ecosystem, finite volume method (FVM) has proven its credibility among other methods such as finite difference (FD) and finite element (FE) method. Lawrence et al. [2] used FVM to determine the air velocity vector within corn and interpolated the nodal velocities in their 3 D FEM model to predict the temperature profiles during aeration.
The FVM uses the conservation principle to represent and evaluate partial differential equations (PDEs) in the form of algebraic equations. The PDEs with the volume integral are converted into surface integrals by the divergence theorem that is evaluated as fluxes on the respective finite volume surfaces. A major advantage of FVM is its adaptability to unstructured meshes that are highly applicable to different shapes and sizes of grain storage silos.
With the advent of computerized platforms, the FVM numerical approach is now widely available in the engineering science discipline of commercial computational fluid dynamics (CFD) software ANSYS FLUENTV R .
The equations that govern the heat, moisture, and momentum transfer are expressed in very general terms and they do not represent the actual stored grain model and the sophisticated aeration parameters. Thus, the CFD codes must be tailored to provide the best fit for the specific aeration model. Thorpe, [3] Ranjbaran et al. [4] and Liu et al. [5] have used similar features to model grain aeration. However, their approach was limited to dimensionality, solution discretization methods and thermo-physical transfer processes that could limit their potential while simulating aeration in a large-scale storage facility.
The model developed by Thorpe [3] was a steadystate 2D model with constant initial and boundary conditions. Besides, the model was not validated by the experimental observations. Ranjbaran et al. [4] developed a 3D transient model with superficial velocity to predict temperature and moisture transfer during the grain drying process. They used secondorder schemes to discretize the governing equations for a total simulation time of 150 min and used a small structure (0.4 m Â 0.4 m Â 0.25 m porous bed) to validate their model. Liu et al. [5] proposed a similar model to predict temperature variation during the grain cooling process in a warehouse (bunker). They assumed grain moisture to be uniform throughout their storage facility to reduce computational complexity. Besides, they used first-order schemes to discretize the governing equation.
The heat and moisture governing equations are linked together with a rate change of moisture in each finite volume that occurs during the aeration process. This shows that the difference in interstitial air and grain moisture could affect the temperature gradient in the stored bulk. The changes in moisture within each finite volume are controlled by the source terms of respective equations that acts as the driving force for the temperature and moisture transfer after the system has attained local equilibrium.
The heat source term considered in the works of Ranjbaran et al. [4] and Liu et al. [5] was defined by the latent heat of vaporization of free water as a function of temperature during the simulation. The heat of vaporization is the energy needed for the phase change process. However, they did not consider the additional energy that must be added to remove the water vapor from the grain surface. This total energy is termed the differential heat of sorption, which is added into the system to liberate the water as vapor from the grain surface. [6] Similarly, released from the system when vapor is adsorbed as water onto the grain surface while aerating with air wetter than the stored grain. [7] This could potentially add to the conservation of thermo-physical balance in each finite volume of the computational domain. Although Cenkowski et al. [8] observed a small difference in latent heat of vaporization and differential heat of sorption at grain moisture above 16.7% wb, adding the sorption parameter would generalize the model's range during aeration drying and or cooling process below the above value. Therefore, it demands the accurate modeling of initial grain conditions in each layer that could contribute to the accurate release or addition of heat of sorption during grain aeration.
Besides, all the previous models have considered superficial velocities to govern the heat and moisture transfer processes in their computational domain. However, superficial velocity is always lower than the actual velocity in the porous medium which could potentially result in slower temperature and negligible moisture changes in comparison to the actual case. [9] Furthermore, the selection of the solution methods to discretize the governing equations can have a significant influence on the accuracy of the simulation results and its associated computational expense.
Therefore, the study aimed to develop a generalized 3D transient heat and moisture transfer model (maintaining optimum accuracy and with the nominal computational time) using ANSYS Fluent that can be used to predict temperature and moisture transfer during the grain aeration process in on-farm silos.
The following objectives were carried out to achieve the above aim: different solution discretization schemes such as pressure-velocity coupling schemes, gradient schemes, spatial discretization schemes, and transient formulation schemes on the prediction error between experimental and model output.
Results were compared with Ranjbaran et al. [4] and Liu et al.'s [5] error profiles to emphasize the process incorporation.
(e) To analyze the variations in differential heat of sorption with respect to the latent heat of vaporization along the stored bulk height.

Experimental procedures
2.1.1. On-farm stored grain facility and data collection The on-farm grain storage facility opted for this study was a flat-bottomed corrugated steel silo located at Balaklava, South Australia, Australia (34 7 0 26.22 00 S, 138 27 0 41.58 00 E) described in the previous work. [10] The on-farm 1000 t silo 12.8 m in diameter and 10.06 m in eave height was filled with 970 t of 2-row feed barley (Figure 1a). The grain level was measured to be 9.91 m high from the silo floor with a 21 angle of repose. The silo floor had a V-shaped aeration floor of 10.4 m in length and 0.9 m in width ( Figure 1c). Detailed information on the aeration ducting configuration is mentioned in the previous work. [9] Four OPI Blue monitoring cables (OPIsystems Inc, Calgary, AB) with Sensirion sensors (error of ±1.8% between 10 to 90% RH with ±0.3 C in temperature) denoted as M1/T1 ('M' for moisture and 'T' for temperature), T2, T3, and T4 cables were installed at the center, north, south-east, and west of the silo as shown in Figure 1d. All the cables except M1/T1 were hanged 2 m away from the silo wall. There was a total of 33 temperature and 9 relative humidity sensors on the four cables with 9 sensors on the M1/T1 cable and 8 each on T2, T3 and T4 cables. The RH sensors on the M1/T1 cable used proprietary equilibrium moisture content (EMC) curves to compute the grain moisture content using the intergranular air temperatures and humidities. All the sensors were at 1.2 m spacings except the bottom sensor which was 1.5 m high from the silo floor. Initial grain temperature and EMC before the fan operation is mentioned in Figure 7.
Before the start of the fan, a stainless-steel truck spear (Graintec Scientific; 1.83 m long and 35 mm diameter; 12 openings) was used to take grain samples near the central cable (M1/T1) to measure the moisture content. Three samples each were taken adjacent to the top two sensors (buried in grain) for correlating the moisture values. Approximately, 10 g of the collected barley samples were placed in the oven at 130 ± 1 C for 20 h to determine the grain moisture content. [11] The initial measured moisture content adjacent to the OPI sensors were 9.7 ± 0.1% wb (top) and 10.5 ± 0.1% wb (second to top), respectively. While the sensors showed 9.9 and 10.8% wb, respectively. Due to the non-feasibility of obtaining actual grain moisture (conventional method) in the large storage silo on an hourly basis, the OPI sensor data were considered for validation processes because of their close resemblance.

Grain cooling operation
A continuous aeration cooling period of 24 h was taken from the previous study. [10] Hourly grain temperature and EMC changes during cooling were taken from the OPI sensors for the field validation of the model. Ambient air conditions during the cooling process were taken from OPI's plenum sensor positioned in the transition duct and are mentioned in Figure 4. The reason behind opting this system was to avoid further calculations on fan warming properties for accurate inlet data. Throughout the study, plenum data will be referred to as the ambient data to provide generality to the situation. To develop a model, it is crucial to adopt certain assumptions to simplify the complex process. Thus, the following assumptions were considered during the aeration simulation: (a) Grain and the inter-granular air in each element are in local thermal equilibrium with each other as described by Thorpe and Whitaker; [12] (b) There is negligible grain shrinkage or expansion in an element taking place during aeration.
According to these assumptions, the transient heat, mass, and momentum balance for the grain and air (combined) were established in each control volume and the following transient-based governing equations were obtained.
2.2.1.2. Momentum transfer model. The momentum transfer model used in this study was taken from the previous study. [9] The governing equations were solved for the physical velocity formulation with incompressible laminar flow in the porous media with two additional resistance source terms to account for the pressure gradient in the domain.
Conservation of mass: Conservation of momentum: where e is grain porosity (decimal), q a is the density of air (kg m À3 ), t is time (s), ! velocity vector (m s À1 ), rp pressure gradient (Pa m À1 ), l is the viscosity of air (kg m À1 s À1 ), (1/a) is the viscous resistance term (m À2 ), C 2 is the inertial resistance term (m À1 ), ! is velocity magnitude (m s À1 ).
The first term represents the pressure gradient in the porous element. The second term represents the diffusive momentum transfer due to the viscous effect. The third term represents the viscous and inertial drag force imposed by the porous walls on the fluid medium and is described by Ergun's equation (Equation (3)).
where DP is the static pressure drop (Pa), L is the element size (m), d p is the grain geometric mean diameter (m). The 1/a and C 2 terms in Equation (2) were expressed considering the viscous and inertial term of Ergun's equation as: viscous resistance component inertial resistance component The anisotropic porous medium components such as variations in porosity and resistance losses were taken from the previous study. [9] The grain porosity was linearly interpolated from 0.3907 at the core to 0.4088 near the periphery. The horizontal resistances were predicted to be 50% of the vertical airflow resistances.
2.2.1.3. Moisture transfer model. The moisture transfer between the grain and the intergranular air was modeled using the User-Defined Scalar (UDS) function with absolute humidity (x) of air as the scalar quantity: where x is absolute humidity (kg water/kg dry air) of air; D eff is the effective diffusion coefficient of moisture (m 2 s À1 ); S m g is the moisture source term due to moisture adsorption in or desorption from grains (kg m À3 s À1 ).
The transient moisture source term (S m g ) during grain cooling was based on Thorpe, (2008) [3] formulation: where @MC @t where K and N are the grain constants and depend upon the grain type (ASABE standards); t is flow time (s); MC is initial grain moisture content (% d.b.); EMC is equilibrium moisture content of grain (% d.b.). Sorption isotherm proposed and validated by Basunia and Abe [13] for barley grain was adopted to determine the grain moisture content: inversion of the above equation gives: where A, B, and C are grain constants; RH is relative humidity of intergranular air (%). The thin-layer constants for barley grain was taken from Basunia and Abe [13] valid at a wide range of temperature and RH: The RH of the intergranular air can be defined as: where p is the vapor pressure of water; p sat is the vapor pressure of water in saturated moist air at any temperature defined as: The value of RH and p sat obtained from Equations (10) and (14), respectively can be used in Equation (13) to obtain the vapor pressure of the intergranular air followed by determining the absolute humidity of the air: where q g is grain's bulk density (kg m À3 ); Ε a and Ε g are air and grain enthalpy (J kg À1 ); K eff is effective thermal conductivity (W m À1 C À1 ); T is the local temperature ( C) and S e g is the energy source term (W m À3 ).
The enthalpy of interstitial air and grain is derived as follows: where C pa (J kg À1 C À1 ) and C pg (J kg À1 C À1 ) are the air and grain's specific heat capacity. The effective thermal conductivity in the porous medium was determined by considering the volume average of the grain and air's conductivity, respectively as follows: where K a and K g are the air and grain thermal conductivity (W m À1 C À1 ). The transient energy source term (S e g ) for grain aeration cooling purpose was adapted from Thorpe, (2008) [3] formulation: where h s is the differential heat of sorption liberated or extracted during moisture adsorption or desorption process (kJ kg À1 ); q t is the true density of grain (kg m À3 ), @MC @t is the rate of grain moisture adsorption or desorption (Equation (8)).
The differential heat of sorption (h s ) was derived from the combination of Clausius-Clapeyron equation and ideal gas law with the consideration of similar grain and free water temperature (based on assumed local thermodynamic equilibrium condition). [14] Then a mathematical expression for the ratio of heat of sorption (h s ) and latent heat of vaporization of free water (h v ) (kJ kg À1 ) considering sorption isotherm (Equation (10)) of grain and saturated vapor pressure of water (Equation (14)) was formulated as described below: h v ¼ 2502535:26 À 2385:76T ð Þ ; ð0 C < T < 65:5 CÞ (22) where p sat is the saturation vapor pressure (Pa); ERH is equilibrium relative humidity (% RH).

Thermophysical property equation of grain
The initial and transient changes in barley grain properties were calculated from the temperature and moisture-dependent formulations as shown below:

Thermophysical property equation of ambient air
The temperature and RH-dependent property equations for humid air 0 C < T < 100 C; 0% < RH < 100% ð Þ for ambient aeration cooling were estimated by the following correlation (Coefficient of Determination: 0.9999): [17] q a ¼ 1:

Stored grain domain
The computational domain was taken from our previous study. [9] A mesh independent study was conducted using 0.25-0.5 m mesh sizes considering static pressure variation and computational time to reach the steady-state convergence. Results showed that a numerical variation of 1% was predicted between 0.35 and 0.25 m sizes while about 6-times higher time was taken by the later size. Therefore, a linear tetrahedral element with 0.35 m mesh size was considered which resulted in 52,788 nodes and 2,93,485 elements. The curvature and proximity mesh sizing and outlet face sizing criteria were considered to generate smaller size elements near the V-shaped aeration floor and at the domain outlet, respectively to improve the mesh quality as mentioned in Figure 2.

Defining grain initial condition
The grain domain was partitioned into 33 zones to account for the temperature and 9 zones for the moisture content, respectively as per the OPI Blue monitoring sensor placements (assigned UDF is shown in Appendix A). An assumption was made that a particular sensor temperature is uniform within a particular zone as specified in Figure 3. T1 zone was assumed to range in a circular radius of 3.2 m (half the radius of the silo). A mathematical discrete function "atan2(Y, X)" was used to define T2, T3, and T4 zones. While moisture (M1) was assumed to vary uniformly across the circular radius of 6.4 m due to a single sensor cable at the bulk center. The height of each zone was constructed such that each sensor placement in the actual silo was equidistant from both sides. The initial grain temperature and moisture values were defined by array of variable. They were denoted as T1 [9], T2 [8], T2 [8] T3 [8], and M1 [9] where the integer in the bracket defined the number of zones along the grain height. All these variables were written under the DEFINE_INIT macro to initialize the moisture content of the grain domain and stored as a user-defined memory variable (C_UDMI). Absolute humidity (xÞ of the domain was determined using Equations (9), (10), (13)- (15) subsequently that was stored as a user-defined scalar variable (C_UDSI).

Defining grain boundary conditions
A volumetric airflow rate of 0.954 m 3 s À1 was observed from a 3-kW aeration fan attached to the on-farm silo. The flow rate and the duct surface area were used to determine the inlet velocity to be 0.051 m s À1 . The no-slip condition was assumed between the wall and grain surface along the silo circumference. Ambient temperature and absolute humidity values were used as the inlet prevailing during the 24 h aeration process (Figure 4). The absolute humidity of the inlet air was determined using Equations (9), (10), (13)-(15) subsequently.

Time-step sensitivity test
The accuracy of the transient formulation depends on the time step size. [19] Though the heat and moisture transfer process in the stored grain ecosystem is a slow process, the too-high time-step value could lead to unphysical results as it directly influences the numerical iterative process of the solver. However, too little could lead to excessive computational time. Based on these considerations and the current domain size, time step sizes of 1, 2, 4, 6, 8, 10, and 12 s were tested for the aeration cooling simulation. Courant number was considered constant due to the assumed incompressible fluid and implicit method usage.     Node-based and Least Squares Cell-based interpolations were examined for Gradient discretization. PRESTO! (PREssure Staggering Option) and secondorder upwind scheme-based interpolation method was used for pressure and momentum discretization, respectively because of the presence of a porous medium. [9,20] While First and Second-Order Upwind schemes for Energy and UDS discretization were considered consecutively for examination.

Transient formulation schemes
In this simulation, both First and Second-Order Implicit transient formulations were examined to analyze the computational time and accuracy in the numerical solution.

Convergence
The absolute criteria for continuity, x, y, and z-velocity and UDS were set at 10E-04 while 10E-06 was set for energy. For making a proper decision on solution convergence, the value of the velocity, energy, and UDS imbalances was monitored at each time step for every iteration (imbalances should be less than 1%) and the volume flow rate at the inlet and outlet boundary section were ensured to stabilize.

Computational procedure
FLUENT's parallel solver (Appendix A) was preferred over serial processing as an advantage for computing a solution to large models by using multiple processes on the same machine. For this reason, the UDFs developed for the aeration model were incorporated with the following compiler directive: "#if PARALLEL (distributes the computing procedure over the number of cores) / Ã Desired Function Ã / #endif (ends the process to generate one single value)." Besides, to avoid the segmentation fault, host processes were excluded by using only #if RP_NODE directive instead of #if PARALLEL. Inlet temperature and absolute humidity values were incorporated in a '.Prof' file for continuous input during the simulation. The stepwise process for model implementation is shown in Figure 5. All the simulations were run in parallel processing with 8 core dual Intel Xeon Gold 5215, 2.5 GHz, 192 GB memory with a bus speed of 2 Â 10:4 GT=s:

Model evaluation and error analysis
After the model was developed and simulated with different parameters and solution methods it was evaluated with real-time temperature (33 data points) and moisture content (9 data points) data taken from the OPI Blue monitoring sensors over each hour in the 24 h operation. The Standard Errors of Prediction (SEP) and Mean Relative Deviation (MRD) between the observed and predicted values were calculated using Equations (31) and (32). The analysis procedure was as follows: first SEP and MRD were calculated for a single sensor point (over the transient basis) followed by averaging with the rest of the sensor points on a single cable. Then, the respective SEP and MRD for each of the cables were averaged (4 cables) and denoted as a single data. Significance (P < 0.05) with respect to time-step and solutions methods as an independent variable was analyzed using Regression method and Single-factor ANOVA, respectively.
where u is the parameter of interest.

Effect of time-step variation
For time-step independence, the total clock time (simulation time), SEP, and MRD results were analyzed for seven different step sizes as shown in Table 1. It is normal that with decreasing time step value the clock time should increase however, results showed that this trend depends upon a specific model and the rate of solution stabilization. All the simulations achieved their first convergence at 162 iterations although they had different step sizes. However, total iterations (Table 1) taken for the simulation showed that at different step sizes, the solution stabilization process is significantly different, eventually affecting the clock time. This could be due to the element size in the domain as the flow properties between cell nodes try to stabilize considering the transient time allotted to the process. Though they were highly dependent, no correlation was observed between them.
In the case of SEP and MRD values for temperatures and MC, significant decreasing error values were observed with decreasing time step sizes except for 1, 2, 4, and 6 s cases with the lowest value for 1 s step size. This could mean that the thermo-physical behavior during aeration cooling was well captured at 1 s in comparison to the 6 s step size during the transient simulation process. Yet, achieving an accurate solution with a nominal computational time can also be an added aspect of a generalized model. Thus, considering all the parameters from the analysis, the 4 s step was selected to keep the calculation cycle. Liu et al. [5] came up with a higher (6 s) time step size for their transient simulation. This could be due to a different computational domain (warehouse) considered for the analysis process.

Effect of solution methods
Simulation results showed that for the 970-t stored barley with 52,788 nodes and 2,93,485 elements, SEP for temperature ranged between 1.31 and 1.42 C among different solution methods. While MRD for temperature ranged between 3.7% and 4.3% (Table 2). Negligible moisture changes (Figure 7) were observed and predicted during the 24 h aeration period justifying its slow cooling front.
The SEP values for temperature and MC obtained with the current model were observed to be similar to that of Ranjbaran et al. [4] while MRD values were observed to be 6.3% (temperature) and 3.8% (MC) less than their reported values. Besides, in all the previous studies, the MRD values between the predicted and experimental data ranged between 3 to 9% for temperature and moisture content. This improvement could be attributed to the inclusion of several factors such as modeling physical velocity, higher-order thermo-physical property equations of air and barley grain, and the inclusion of differential heat of sorption.
Besides, the deviations observed in temperature and moisture values could be due to the following reasons: a. Assumption of uniform barley grain moisture (across the circular radius of 6.4 m) and temperature (near the peripheral region) across the thin layer that could have altered the release of the heat of sorption during the aeration cooling process. Figure 6 explains the potential difference in temperature that could arise between the assumed and the actual layer. The initial barley grain conditions for each sensor zone were modeled assuming equal height from the silo floor (for sensors at 1.5 m from the floor) and subsequent layer along the cables. However, during the loading process, grain forms a truncated cone ( Figure 6) above the floor after a typical day of filling. This could potentially result in different temperature and or moisture values in the peripheral zones.
Thin-layer constants (K and N) may not be so accurate for the 2-row feed barley grain variety.
Though the aeration model showed reasonably low error, there is a need to determine the solution  methods that would be useful to run a long-term aeration simulation for a typical storage period. Table 2 shows a complete overview of the solution method analysis. It was observed that different solution methods had taken noticeably different computational times (P < 0.05) except between two gradient methods with SIMPLE algorithm to carry out an actual aeration process of 24 h. Besides, P > 0.05 was observed for temperature within the pressure-velocity and gradient schemes however, P < 0.05 was detected for spatial and transient formulations.
Ranjbaran et al. [4] used 0.01 s to simulate the heat and moisture transfer for 150 min. They used a second-order scheme to evaluate the transient part of the discretized equation. This implementation could be justified as they used a small simulation time of 150 min for a pilot-scale bin. While with similar solution schemes (SIMPLE, Least Square Cell-based and second-order schemes for spatial discretization) as that of the above study (highlighted in bold in Table   2), the current aeration model took almost twice (% 190% more) the simulation time.
Besides, Ranjbaran et al. [4] used a second-order scheme for pressure discretization while Liu et al. [5] used a first-order scheme for momentum and pressure discretization that is not recommended due to discontinuous pressure gradients imposed by the porous medium in each element. These schemes could worsen the solution stabilization leading to slower convergence that is taken care of by the PRESTO! scheme due to its staggering continuity balance about the finite volume face. [21] This has been further validated in this study as shown in the simulation time data.
Though ANSYS Inc. [22] suggests using the COUPLED algorithm for the pressure-velocity scheme when large time steps are used, the transient stored grain aeration analysis showed similar temperature and moisture error as compared to SIMPLE (Table 2). However, the former scheme took almost twice (% 190% more) the time (P < 0.05) with First-order spatial and transient schemes, respectively. This difference could be attributed to the presence of a porous medium in the domain.
The discontinuous pressure gradient enforced by the porous domain could execute convergence issues (within momentum and pressure-corrected continuity equation) with non-linear COUPLED algorithms thereby taking more time (P < 0.05) to converge than that of the SIMPLE algorithm.
Similar results were observed between the two gradient schemes with the Least-Square Cell-based method taking lesser time than the Green Gauss Node-based method yet they were insignificant (P > 0.05). This could be due to the complex process of evaluating the face value of the respective scalar quantity (temperature and or absolute humidity) in the finite volume. In the former scheme, the solution at the finite volume face is assumed to vary linearly from the cell centroid however, the solution is computed by the arithmetic average of the nodal values of each face of the finite volume in the later scheme. [22] The exploration of Green Gauss Node-based gradient could be beneficial for larger storage silos (>1000 t) as their mesh size would be >0.35 m as opted in this study. This was further validated by a previous study [23] where they observed a 75% increase in the mesh size when transitioned from 100 to 1000 t storage silo. Furthermore, there could also be the influence of the shape as the 100t silo had a hopper-bottomed section. Thus, considering a nominal computational time, SEP and MRD values from the simulations, the optimum solution methods confirmed for the Pressure-velocity scheme, Gradient, Energy discretization, UDS discretization, and Transient formulation were SIMPLE, Green Gauss Node-Based, First-order, First-order, and First-order schemes, respectively. Therefore, from the results presented in Table 2, it can be confirmed that the laminar flow in a tetrahedral element of grain porous medium aligns reasonably with the mesh size. Justifying the acceptability of first-order upwind discretization for energy and UDS that resulted in insignificant numerical discretization error (numerical diffusion) at each iterative process. Figure 7 shows the thermo-physical phenomenon with optimized solution methods. The model closely predicted the flow of the warm and cool front from the subsequent layers. Among all the four cables, predicted data along the T1 cable sensors were more scattered when compared with the actual data. Due to the positioning of the T1 cable along the core of the stored bulk, grain at 1.5 m was observed to undergo a slower cooling rate than its predicted counterpart and subsequently, its effect was forced upwards. Barley grains at 11.2 m high (T1 cable) and 9.9 m high (T2, T3, and T4 cables) exhibited the highest increase in temperature due to the incoming warm front from the grain below the top layer. A temperature lag was observed at 9.9 m high on the T1 cable whose effect was transferred to the grain at 11.2 m from the floor. This could be due to lower airflow justifying the random barley grain property distribution in the peaked grain volume. Ambient RH was predicted to have a reverse effect on the barley grain moisture (at 1.5 m high from the floor) due to lower absolute humidity than the grain (Figure 8). Similarly, subsequent upper layers were affected by humidity front resulting in non-uniform adsorption and desorption process in the on-farm silo. The above simultaneous process affected the heat of sorption liberated during the aeration process. Figure 9 showed that the differential heat of sorption was higher than the latent heat of vaporization of pure water at any depth of the grain bulk. The difference ranged between 25 and 45 kJ kg À1 of water adsorbed by the grains with the lower value at 1.5 m and upper values at 11.2 m from the silo floor. This justified that the energy that holds the water together in the liquid phase is lower than the energy of binding between the sorption (hydrophilic) sites and the water molecules. Besides, studies on the grain molecular level have stated that at low moisture, sorption sites are highly available where adsorption takes place at the monomolecular layer of the grain kernel while at high moisture, the strength of water-binding decreases. [24] McMinn and Magee [25] suggest that this could be due to the changes in the geometry of polymers during sorption that ultimately affects the activation of sorption sites on the grain particles. These molecularlevel changes have been attributed to closely resemble the actual temperature changes observed in the subsequent upper layers (Figure 7) in the stored bulk. Figure 10 shows the effect of ambient conditions on barley grain within 1.5 m high from the inlet. Grain moisture increased up to 13.4% (wb) at the interface of grain and inlet trench thus reaching equilibrium with the ambient air. However, the degree of equilibrium state gradually decreased up to 0.04 m. This phenomenon is per the psychrometrics of air yet gets affected by the grain condition as shown at 1.5 m high and above.

Aeration cooling effects with optimized methods
Quantitative data on the distribution of wetting front near air entry in on-farm silos are scarce, however, with the help of computational modeling stored grain managers can be provided with related information. This will eventually contribute to managing the aeration control system for preserving grain moisture within the acceptable limit and if possible, contribute to the return value of dry grains.
Besides, in this paper, an analysis conducted using the current model's output with previous research  showed that even with a similar numerical method, the solution discretization schemes highly differ when a model's efficiency is scaled up to an actual scenario. This further sets-up the stage to scrutinize the effect of scheme switching (First-to-Second order blending) between iterations to see how the convergence affects the solution stability.
Such a model could be a reliable foundation of information for optimizing fan control strategies for aeration cooling management practices.

Conclusion
A 3 D transient model to predict temperature and moisture transfer during on-farm stored grain aeration cooling was successfully developed using the finite volume-based ANSYS Fluent tool. Compiler directives were successfully incorporated to compute solutions for large-scale models. The incorporation of the differential heat of sorption and physical velocity formulation lessened the statistical error between the predicted and the observed scalar components of the aeration model in comparison to previous models. The following are concluded from the above simulation: A 4 s time-step size was optimized for the aeration cooling process.
a. Temperature error profile showed significant fluctuations with spatial and transient formulations while no significance was observed within the pressure-velocity and gradient schemes. Similarly, simulation time showed no significance between Least Square Cell-and Green Gauss Nodebased schemes. b. COUPLED algorithm was observed to take almost 190% more time than SIMPLE algorithm. While time taken with Second-order simulation were marginally higher than its First-order counterpart. Thus, it was concluded that opting a Secondorder and COUPLED schemes are not efficient for aeration simulation in large storage silos. c. The following solution methods were optimized for aeration cooling simulation considering the minimum computational time and high accuracy for Pressure-velocity scheme: SIMPLE; Gradient: Green Gauss Node-Based; Energy discretization: First-order; UDS discretization: First-order; and Transient formulation: First-order. SEP for temperature and moisture ranged between 1.1-1.7 C and 0.3% wb, respectively. While MRD for temperature and moisture ranged between 3.1-4.4% and 2.2%, respectively. d. The predicted differential heat of sorption was found to be higher by 25-45 kJ kg À1 than the latent heat of vaporization during the aeration cooling process.
Such a model could be a reliable foundation of information for optimizing fan control strategies for aeration cooling management practices. Future work by the authors will look to integrate a fan (inlet) control system using user-defined functions in CFD that can help to optimize aeration strategies for any location.