A physical model to predict climate dynamics in ventilated bulk-storage of agricultural produce

International Journal of Refrigeration 30 (2007) 195e204 www.elsevier.com/locate/ijrefrig

A physical model to predict climate dynamics in ventilated bulk-storage of agricultural produce L.J.S. Lukasse*,1, J.E. de Kramer-Cuppen, A.J. van der Voort Wageningen University and Research, P.O. Box 17, 6700 AA, Wageningen, The Netherlands Received 1 April 2005; received in revised form 30 January 2006; accepted 13 March 2006 Available online 5 June 2006

Abstract This paper presents a physical model for predicting climate dynamics in ventilated bulk-storage of agricultural produce. A well-ordered model presentation was obtained by combining an object-oriented zonal decomposition with a process-oriented decomposition through matrixevector notation. The objective of this paper is twofold: (1) to present the modelling procedure and (2) to present the resulting model, validated with real data. The model is a suitable simulator to assess potential effects of changes in ambient climate, design, and controller tunings. The model predictions fit well to extensive real data from three different cases. The good fit for all three cases was achieved with the five calibration parameters calibrated on the basis of data from one case only.  2006 Elsevier Ltd and IIR. All rights reserved. Keywords: Fruit; Vegetable; Storage; Ventilation; Modelling; Variation; Temperature; Humidity

Mode`le physique permettant la pre´vision de la dynamique climatique dans l’entreposage des fruits et le´gumes en vrac Mots cle´s : Fruit; Le´gume ; Stockage ; Ventilation ; Mode´lisation ; Variation ; Tempe´rature ; Humidite´

1. Introduction Many models to predict climate conditions in agricultural storage facilities are available in the literature ([1,3e10]). Unfortunately in many papers presenting models of complex systems, like climate in storage facilities, the complex notation hampers the understanding of * Corresponding author. Tel.: þ31 317 475301; fax: þ31 317 475347. E-mail address: [email protected] (L.J.S. Lukasse). 1 IIR Member of Commission D2. 0140-7007/$35.00  2006 Elsevier Ltd and IIR. All rights reserved. doi:10.1016/j.ijrefrig.2006.03.011

interactions between the different model states. This paper presents a physical model to predict climate dynamics in ventilated bulk-storage of agricultural produce. This paper enriches the available literature on modelling distributed climate dynamics due to the model’s completeness, the well-ordered model presentation, and the model validation with real data. The ordering was achieved by combining an object-oriented zonal decomposition according to [1] with a process-oriented decomposition through matrixe vector notation as presented in [2]. The objective of this paper is twofold: (1) to present the modelling procedure and (2) to present the resulting model,

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

196

Nomenclature Style bold BOLD italic

non-capital bold symbols represent a vector bold capital symbols represent matrices italic symbols represent a scalar

Subscripts b bulk a component air pr component produce h concerning heat CO2 concerning carbon dioxide w concerning water Variables h dynamic viscosity of air (Pa s [1.71  105]) 3 porosity, cq. void fraction ((m3 air) (m3 bulk) [0.34]) f air flow volume ((m3 air) h1) l heat conduction coeff. of potato (W m1  C1 [0.55]) Fb diagonal matrix with dilution rates for the fluent part of x (h1) ra air density at 5  C (kg m3 [1.27]) rbulk bulk density (kg (m3 bulk) [670 for potatoes]) rpr produce density (kg (m3 produce) [1014 for potatoes]) Aa,pr specific surface area of produce ((m2 produce) (m3 bulk) [60 for average potatoes, Eq. (12)]) Ab,wall wall area in contact with bulk (m2) aw produce water activity ((-) [1]) cp,a specific heat content of dry air (J kg1  C1 [1000]) cp,pr specific heat content of produce (J kg1  C1 [3600 for potatoes]) d diameter of (spherical) produce (m) Er oxidation energy of glucose (ATP þ heat) (J (mole glucose)1 [2.816  106]) K ˛Rnm , gain matrix relating rðxÞ˛Rm to x˛Rn (-) Ka,pr water vapour transfer coefficient ((g water) m2 Pa1 h1 [8.4  104]) Kevap evaporation coeff. ((kg air) (m2 produce) h1 [0.14]) Kwall heat transfer coefficient of wall (W m2  C1 [0.5])

validated with real data. The model predicts temperature, relative humidity, CO2 and weight loss in bulk-storage using input information on ambient climate and design, operation and produce parameters. Bulk-storage is common for potatoes, onions, grain and rice. Moreover in storage of fruit or flower bulbs in bins with a forced airflow through the bins,

ls M min(x,y) pdry air pa p1,resp p2,resp p3,resp qevap Ra r(x) rb,h,apr rb,resp rb,evap rb,wall Rpr RQ Ta Text Tskin v Vb ws Xa,sat xib xb;CO2 ;a xb,h,a xb,h,pr xb,m,pr xb,w,a xin

store length (m) molar mass (g mole1) minimum of x and y pressure of dry air (Pa [105]) partial vapour pressure in air (Pa) respiration parameter (Eq. (17)) (W tonne1 [10.75]) respiration parameter (Eq. (17)) (W tonne1  C1 [1.56]) respiration parameter (Eq. (17)) (W tonne1  C2 [0.17]) latent heat of evaporation (J g1 [2490]) resistance to heat transfer from produce to air (m2  C W1) ˛Rm , vector of process rates (-) rate of heat transfer from air to produce (J (m3 air)1 h1) respiration rate (J (kg produce)1 h1) rate of moisture evaporation from produce ((g water) (kg produce)1 h1) heat exchange through walls (J (m3 air)1 h1) heat flow resistance within a single produce (m2  C W1) respiration quotient (mole CO2 produced per mole O2 consumed [1]) air temperature (  C) external (ambient) temperature (  C) produce skin temperature (  C) air velocity (m s1) volume of bulk layer i (m3 bulk) store length (m) saturated vapour concentration ((g water) (kg air)1) ˛Rn , state vector x of zone i in bulk (b) CO2 concentration in the bulk (%) heat content of air in the bulk (J (m3 air)1) heat content of produce in the bulk (J (kg produce)1) cumulative produce mass loss in the bulk (g (kg produce)1) water vapour content in the bulk ((g water) (kg air)1) state vector x of influent air

the local climate conditions within the bins exhibit very similar characteristics, albeit on a smaller scale. The model was developed within a joint research project. In that project the model was used to analyze the potential effects of control algorithms for automatically anticipating the weather forecasts. Yet the application potential ranges

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

further: it can serve as a simulator to assess potential effects of changes in ambient climate, design, and controller tunings. The outline of this paper is as follows. Section 2 explains the modelling procedure. In Sections 2.1e2.6 the procedure is applied to a bulk-storage facility for potatoes. In Section 3 the resulting model is simulated and calibrated to data collected from real storage facilities. Finally in Section 4 this paper’s conclusions are drawn.

2. Modelling procedure The modelling procedure consisted of six steps: 1. Analysis of flow characteristics and subsequently the selection of zones within which spatial distribution is of minor importance. 2. Identification of relevant model states x (Eq. (2)). This obviously relates to the purpose of the model. For example, when the model is meant to study temperature gradients in regular atmosphere stores then one may well omit O2 from the model state, while if one intends to study the O2 pull-down speeds in Controlled Atmosphere stores then O2 must be in the model state vector x. 3. Identification of relevant processes r(x) (Eq. (2)), i.e. to define the elements of column vector r in Eq. (2). It depends on the modelling purpose which processes are relevant. For example, if one wants to study weight loss issues one will have to model the evaporation process revap, while if one is only interested in steady state temperature-distribution in high-humidity stores then the effect of evaporation may well be neglected. An analysis of the sensitivity of model states to the candidate processes may help to identify the most relevant processes. 4. Definition of the coefficients in gain matrix K (Eq. (2)). It is recommended to adopt a natural sign-convention: if element K(i,j) < 0 then the i-th state diminishes due to process j, i.e. consumption, while K(i,j) > 0 indicates accumulation of state i due to process j, i.e. production. 5. Definition of explicit algebraic output equations y ¼ g(x) mapping model state x to model output y, e.g. mapping absolute humidity (g kg1) to relative humidity (%) or heat content (J kg1) to temperature (  C). 6. Mutual connection of all zones involved. Within one zone in steps 2e4 the system dynamics were modelled using a process-oriented decomposition through matrixevector notation as presented in [2]. This notation is the best available way to present models of complex dynamic zones with perfect mixing, because it facilitates the reader tracing the interactions between multiple processes and states. The heat and mass balance for a state x within

197

a defined zone boundary can be separated in transport and reaction terms: accumulation ¼ input  output þ reactions

ð1Þ

The transport terms input and output depend on physical system characteristics. When modelling a specific zone one first needs to define the zone boundary, and define the transport terms through the zone boundary as input and output. For ventilated systems, like storage facilities, the dominant transport term for a specific perfectly mixed zone is: inpute output ¼ f(xin  x), in which the air flow f is generated by fans or natural convection. Eq. (1) in matrix-notation: x_ ¼ fðxin  xÞ þ K  rðxÞ

ð2Þ

In Eq. (2) the matrixevector product K  r(x) comprises the term reactions in Eq. (1). For systems with spatial gradients in x, like bulk-stores, Eq. (2) is an approximation as the assumption of prefect mixing is violated. The larger the zones the larger the approximation error. Selection of the number of zones is done by first conducting a timeconsuming simulation of a step response for an extremely large number of zones, and subsequently simulating for a reduced number of model zones and balancing the reduction in computation time against the goodness of fit to the original simulation. As with any model, upon completion, the model output must be compared to real data. This comparison may necessitate the modeller to go back to one of the preceding steps in order to define the input variables so that the model outputs better represent the system under investigation. As this is generic for any modelling application, it is not regarded as a characteristic step in this paper’s modelling procedure.

2.1. Step 1: spatial zoning of a bulk-storage facility Fig. 1 depicts the common outline of potato storage facilities. Fans take in a mixture of ambient air and return air and blow the mixture into the air lock. The mixing ratio is controlled by manipulating the position of the hatches. Through a grated floor or air ducts the air enters the bulk at the bottom and flows vertically towards the head space, meanwhile exchanging heat, moisture and gases (of which CO2 is most relevant) with the bulk. Part of the air is discharged from the head space to the ambient, while the remainder is drawn into the fans again. For more details on potato stores see e.g. [11]. Temperature gradients within the storage facility primarily occur vertically inside the potato bulk during cooling down or warming up, due to the flow-through nature of the bulk in combination with the huge difference in heat capacity between the fluid (air) and solid (produce). To facilitate simulation of these one-dimensional in-bulk temperature gradients the produce bulk was modelled as a series connection of multiple horizontal layers (zones). The remaining

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

198

external climate Head space

hatches evaporator fans

Bulk

air Lock Fig. 1. Common outline of potato storage facilities.

zones were designed in analogy with existing physical boundaries within the store (Fig. 2). For reasons of space limitations only the model for one zone ‘one bulk layer i’ is presented in detail in the next sections. The model equations for the zones ‘air lock’ and ‘head space’ can be deduced from the bulk layer model primarily by omitting the states ‘produce temperature’ and ‘produce mass’ as well as the processes rb,h,apr, rb,resp and rb,evap (Eq. (6)). 2.2. Step 2: applied to a bulk layer: selection of model states It was felt that the current situation, state, in a specific location in the bulk, is sufficiently known when produce heat content (or temperature), air heat content, humidity, CO2 and cumulative weight loss were known. The direct reason to include heat content, humidity and CO2 in the model state was that this model’s prime application was to evaluate the effect of control strategies. The current generation of climate

controllers for potato stores actively manipulates temperature, humidity and CO2. Hence to allow evaluation of simulated control strategies the effect on temperature, humidity and CO2 had to be predicted. The reason to discriminate between air and produce temperature was the large difference in heat capacity and the consequential large temperature differences in non-steady conditions. Cumulative weight loss is not a climate characteristic, but an easy-to-model produce quality characteristic affected by the climate control strategy. Therefore this produce quality characteristic was included in the model state as well. As a result a bulk layer (Fig. 2) was defined as a zone described by five states. Three of them described the fluent air while the remaining two concerned the solid potato. The state xib of bulk layer i is: 2 i 3 xb;h;a ½J ðm3 airÞ1  i 6 xb;w;a 7 ½ðg waterÞ ðkg airÞ1  6 i 7 7 ½% xib ¼ 6 ð3Þ 2 ;a 7 6 xb;CO i 4 xb;h;pr 5 ½J ðkg produceÞ1  i xb;m;pr ½g ðkg produceÞ1  The states are influenced by transport phenomena (air flow, heat and mass transfer) and a number of processes. In a matrixevector notation this is modelled by i i x_ ib ¼ Fb ðxi1 b  xb Þ þ Kb  rb ðxb Þ

ð4Þ

The transport term Fb is a diagonal matrix, instead of a scalar f, because only the fluent part of the states in a zone is directly affected by the physical transport terms: 2 3 1 0 0 0 0 60 1 0 0 07 7 f6 6 0 0 1 0 0 7 ½h1  Fb ¼ ð5Þ 7 3Vb 6 40 0 0 0 05 0 0 0 0 0 2.3. Step 3: relevant processes r(x)

external climate

Head space

Bulk, layer n

Bulk, layer 1 air Lock (+duct) Fig. 2. Zonal decomposition of a potato store for modelling purposes.

The identification of relevant processes requires physical knowledge of the system involved. One needs to select a priori the physical processes which are regarded relevant to the model states. Of course this selection may be reconsidered if the model is falsified when confronted with real measurement data. A simple steady state heat balance analysis shows that for normally insulated potato stores the two dominant terms in a store’s heat balance are heat exchange due to conduction through the walls rb,wall and autonomous heat production by respiration rb,resp. Of course respiration is part of complex metabolic networks involving thousands of biochemical processes inside the potato tuber, however, this process can be simplified in the model to a lumped process rb,resp. The choice made in Section 2.2 to distinguish between produce and air temperature necessitated the extension of the model with rb,h,apr, convective heat exchange between air and produce. The last process included was evaporation of

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

moisture from the produce rb,evap. It was included in the model because it is the main cause of weight loss and it affects the air humidity (both model states). As a summary of the preceding description of processes the vector of relevant process rates was given by T rb ¼ ½ rb;wall rb;h;apr rb;resp rb;evap 

Rpr

rb;wall ¼ Kwall Ab;wall ðTb;a  Text Þ3600=ð3Vb Þ ½Jðm3 airÞ1 h1  ð7Þ where i xb;h;a ½  C ra cp;a

The value for Kwall, the heat transfer coefficient of the wall, was taken from [12]. 2.3.2. Heat transfer from air to produce (step 3) The process of heat transfer from air to produce, rb,h,apr, is a series connection of two processes: external transport from air to produce skin, and internal transport from skin to deeper inside the produce (Fig. 3). The driving force is the difference between air temperature (Ta) and volume average produce temperature (Tpr) ([13]), while the heat flow resistance is a series connection of the internal resistance Rpr and the external resistance Ra. The overall heat transfer process was modelled as rb;h;apr ¼

  1 Aa;pr Ta  Tpr 3600=rbulk Ra þ Rpr ½J ðkg produceÞ1 h1 

where dc d ½m2  C W1  Rpr ¼ ¼ l 8l

air (Ta)

produce (Tpr) internal transport

The four following sections will each discuss one individual element of rb (step 3). The consequence of these processes for the states, i.e. for the contents of the accompanying column of Kb (step 4) will be discussed in Section 2.4.

Tb;a ¼

Ra

ð6Þ

2.3.1. Heat conduction through walls (step 3) Both the temperature of the produce and the air in the storage facility determine the heat transport through the store walls. For the overall model outcome it made little difference how this mutual distribution was modelled. Therefore the model could be simplified in this perspective and it was assumed that all heat comes from the air in the bulk:

ð8Þ

ð9Þ

199

transfer from air to produce

Fig. 3. Heat transfer from air to produce.

v ¼ f=ðls ws 3600Þ ½m s1 

ð11Þ

2

pd ð1  3Þ ½ðm2 produceÞ ðm3 bulkÞ1  1 3 pd 6 i i xb;h;pr xb;h;a ½  C; Ta ¼ ½  C Tpr ¼ cp;pr ra cp;a

Aa;pr ¼

ð12Þ

ð13Þ

In Eq. (9) dc is the characteristic diameter of a sphere ( ¼ d/8). That is for a spherical object the time constant of a cooling or warming curve with dc ¼ d/8 equals the average time constant of the curves resulting from an exact solution of the partial differential equations solving the exact temperature transients within a spherical object. Note that Eqs. (9) and (12) assume a spherical produce, then A ¼ pd2 and dc ¼ d/8 ([13,14]). Deviations from this spherical shape will cause small prediction errors, for example a cylinder has dc ¼ d/6. The arctan-function in Eq. (10) was taken to facilitate a smooth transition from the natural convection regime at low air velocities (Ra,NC) to the forced convection regime at high air velocities (Ra,FC): Ra;FC ¼ 0:018v0:59 ½m2  C W1 

ð14Þ

0:49

ð15Þ

Ra;NC ¼ 0:024v

2 

½m

1

CW 

Eqs. (14) and (15) were derived from the common dimensionless heat transfer correlations for packed beds. Eq. (14) holds for Re < 180 and 15 for Re > 180, see e.g. [15]. Although Ra is physics-based its exact value is rather uncertain. 2.3.3. Respiration (step 3) Respiration is the oxidation of organic matter by living organisms in order to generate energy. In storage of fruit and vegetables the respiration process consumes organic

     11 180h 0:49 0:024v 0:5 arctan 100 v  þ 0:5 B B CC B   CC ½m2  C W1     ra d Ra ¼ minB @1; @ AA 180h þ 0:5 þ 0:018v0:59 0:5 arctan 100 v  ra d 0 0

ð10Þ

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

200

materials within the produce and produces heat and CO2. The chemical reaction equation is rb;resp

C6 H12 O6 þ 6O2 ! 6CO2 þ 6H2 O þ Er

ð16Þ

where Er ¼ oxidation energy of glucose (ATP þ heat) ¼ 2.816 MJ (mole glucose)1. The temperature dependence of rb,resp was modelled as ([4])     rb;resp Tpr ¼ p1;resp  p2;resp Tpr þ p3;resp Tpr2 103  3600 ½J ðkg prÞ1 h1 

ð17Þ

where Tpr follows from Eq. (13). The temperature dependence of rb,resp(Tpr) for the parameter values listed in the nomenclature corresponds to the data in both [16] and at www.postharvest.com.au. Although temperature is an important factor, there are more factors seriously affecting the rate of respiration. Most of them have to do with the physiological stage of development: e.g. ripeness and wound healing. These factors are hard to measure and vary depending on year, store, season, cultivar, etc., hence causing a wide uncertainty interval for rb,resp (þ or 50% is not uncommon, see e.g. [17]). Eq. (17) assumed that the CO2 concentration does not reach levels at which it inhibits the respiratory activity. In usual potato storage regimes this assumption is satisfied. 2.3.4. Evaporation (step 3) The driving force for evaporation is the vapour concentration difference between the produce skin area and the surrounding air. When the moisture content in the boundary layer of the produce is higher than in the surrounding air, water will evaporate from the produce to the surrounding air. This will usually be the case as respiration causes the produce to be slightly warmer than the surrounding air. Also, the boundary layer inside the skin will always be saturated with water vapour. With aw ¼ 1, the rate of evaporation was given by rb;evap ¼ Kevap

Aa;pr ðXa;sat  xb;w;a Þ½ðg waterÞðkg produceÞ1 h1  rbulk ð18Þ

Here, the driving force is written as a vapour concentration, because this enables us to use a model state. The saturated vapour concentration Xa,sat depends on the temperature of the skin. Like in [1] it was modelled by the Antoine equation for water vapour: pa;sat e Xa;sat ¼ ¼ 162

ð23:5Tskin3991þ234Þ 162

½ðg waterÞ ðkg airÞ1 

ð19Þ

The above expression causes some minor errors in the regions of T > 100  C, but not when 0 < T < 15  C, the

region relevant for storage of perishable produce. The number 162 in Eq. (19) is an approximate conversion from Pa to g/kg derived from the ideal gas law: nv RðT þ 273:15Þ V 1 ¼ Xa ra RðT þ 273:15Þ ½Pa Mw

pa ¼

ð20Þ

At T ¼ 5  C, the usual potato store temperature, ð1=Mw Þra RðT þ 273:15Þ ¼ 162. The evaporation coefficient Kevap in Eq. (18) was derived from the water vapour transfer coefficient Ka,pr [(g water) m2 Pa1 h1)]. Kevap depends on the produce type (esp. peel characteristics) and air velocity. Despite the theoretical dependence on air velocity, reliable quantitative information is not available. Therefore Kevap was modelled as a function of the produce-specific Ka,pr (8.4 (g water) m2 Pa1 h1), based on values in [11] and [18]: Kevap ¼

Ka;pr Mdryair pdryair MH2 O

¼ 0:14 ½ðkg airÞ ðm2 produceÞ h1 

ð21Þ

The produce skin temperature Tskin in Eq. (18) was given by Tskin ¼

Ra Tpr þ Rpr Ta  ½ C Ra þ Rpr

ð22Þ

with Ra and Rpr as defined in Eqs. (9) and (10), and Ta and Tpr given in Eq. (13). Eq. (22) was based on the assumption that the internal and external heat flow in a single produce are always in equilibrium, as in that situation: Tpr  Tskin Tskin  Ta ¼ Rpr Ra

ð23Þ

2.4. Step 4: gain matrix K The complete reactions term Kb  rb (Eq. (4)) was given by Eq. (24). The columns of Kb will be explained beneath. 2 3 0 1 rbulk =3 0 r bulk 60 7 0 0 6 3ra 7 6 7 K b  rb ¼ 6 0 7 0 9:3  103 0 6 7 40 1 1 qevap 5 0 0 2:6  105 1 3 rb;wall 6 rb;h;apr 7 7 6 4 rb;resp 5 rb;evap 2

½J ðm3 airÞ1 h1  ½ðg waterÞ ðkg airÞ1 h1  ½% h1  ½J ðkg produceÞ1 h1  ½g ðkg produceÞ1 h1 

ð24Þ

The 1st column in Kb accounts for the effects of rb,wall, heat leakage through walls, on the model states. In Section

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

2.3.1 it was assumed that heat leaking through the walls affects only the air temperature in the bulk, hence the 1st column of Kb is [1, 0, 0, 0, 0]T. The process of heat transfer from air to produce, rb,h,apr, obviously affects the states i i xb;h;a and xb;h;pr , hence the 2nd column of Kb has two nonzero elements. The term rbulk/3 in Kb(1,2) performs the necessary unit conversion from [J (kg produce)1 h1] to [J (m3 air)1 h1]. The 3rd column of Kb accounts for the effects of rb,resp on the model states. As can be seen from Eq. (16) respiration generates heat, water and CO2 within the produce. The water can only leave the produce by evaporation. CO2 volatilizes to the air, thus increasing xb;CO2 ;a and reducing xb,m,pr. Hence rb,resp affects the 3rd, 4th and 5th model state. The remaining issue was to determine the values of the accompanying gain coefficients Kb(3,3), Kb(4,3) and Kb(5,3). The mass loss resulting from respiration follows from the difference in molar mass between CO2 (mass produced) and O2 (mass consumed) and the rate of respiration: respiratory mass loss ¼

ðRQMCO2  MO2 Þ6 rb;resp Er 5

¼ 2:6  10

The model output equations map the model states x to the model outputs y. In this application they were used to convert the heat contents xb,h,a and xb,h,pr to temperatures Tb,a and Tb,pr, respectively. 2 i 3 Tb;a i 6 xb;w;a 7 6 i 7 7 yib ¼ 6 2 ;a 7 6 xb;CO i 4 Tb;pr 5 i xb;m;pr 3 2 1 0 0 0 0 ½  C 7 6 ra cp;a 7 6 ½ðg waterÞ ðkg airÞ1  7 60 1 0 0 0 7 i 6 x ½% ¼ 60 0 1 0 07 7 b  6 1 ½ C 7 60 0 0 0 4 cp;pr 5 ½g ðkg productÞ1  0 0 0 0 1 ð27Þ

 rb;resp ½g kg

1

h 

hence Kb(5,3) ¼ 2.6  105. The respiratory CO2 production is given by RQMCO2 6rbulk rb;resp Er rCO2 310

¼ 9:3  103  rb;resp ½% h1  ð26Þ 3

2.5. Step 5: output equations

2.6. Step 6: mutual connection of all zones involved 1

ð25Þ

respiratory CO2 production ¼

201

hence Kb(3,3) ¼ 9.3  10 . The unit of rb,resp equals the unit of x_ ib;h;pr and respiration ‘produces’ heat, therefore Kb(4,3) ¼ 1. The 4th column of Kb accounts for the effects of rb,evap on three of the model states: (1) Kb(5,4) e evaporated water is weight loss, (2) Kb(2,4) e evaporated water humidifies the air, and (3) Kb(4,4) e evaporation consumes sensible heat. Kb(5,4) ¼ 1 as the unit of rb,evap equals the unit of x_ ib;m;pr while evaporation ‘consumes’ water from the produce. Kb(2,4) ¼ rbulk/(3 ra) because each gram of evaporated water from the produce ‘produces’ a gram of water vapour in the air, while the units need to be converted from [(g water) (kg produce)1 h1] for rb,evap to [(g water) (kg air)1 i h1] for xb;m;pr . Kb(4,4) ¼ qevap because the quantity of heat ‘consumed’ to evaporate a gram of water equals the latent heat of evaporation (J g1), while no further unit conversion is required. Actually, evaporation draws energy from the produce in two ways: (1) the evaporated moisture cools down from Tpr to Ta and (2) sensible heat is transformed into latent heat qevap by the phase change of water from liquid to gas. This second effect is by far the largest effect, and therefore the first effect is neglected.

The mutual connection of zones comes down to the correct accounting for transport phenomena across zone boundaries. The connection of zones with equal state vectors was straightforward (Eq. (4)). Connecting zones with different state vectors required some precautions in defining the transport terms. For example, the dynamics of the 1st bulk layer were described with: 2 3 1 0 0 60 1 07 7 f6 6 0 0 1 7$xl  Fb $x1 þ Kb $rb ðx1 Þ x_ 1b ¼ b b 6 7 3Vb 4 0 0 05 0 0 0 The above equation differs from Eq. (4) because the definition of xl, the air lock state, differs from the definition of xb: xl ¼ ½ xl;h;a xl;w;a xl;CO2 ;a T . A very similar formulation was applied to the head space zone, because the definition of xh differs from the definition of xb. 3. Model simulation and calibration The model was programmed in a C-mex S-function, within the Matlab/Simulink platform v6.5 (www.mathworks. com). Integration of the stiff set of differential equations was managed by standard Simulink routines. Accurate integration was verified by setting the maximum integration step to an extremely small value, and verifying that the resulting simulation output was virtually equal to the simulation output achieved with the default integration settings (min. and max. step size, abs. and relative tolerance). In order to verify/falsify the model, experimental data were collected from two full-scale bulk-stores with potatoes

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

202

in The Netherlands in the storage season of ’97/’98. For one of the stores, measurements during the last 40 days of the previous season were also available, resulting in three datasets. A model in which the bulk was lumped in one-layer was fitted to the data. The choice to fit a one-layer model was based on the objective to use the model as an internal model for a non-linear Model Predictive Controller, in which highdimensional or complex models are disadvantageous. In a first comparison the difference between model output and measurements was pretty large, as can been seen in the left top graph of Fig. 4. In Fig. 4 (top left) the simulated produce temperature (dashed) is shown for the original parameters. The measurement series (solid) actually is the mean of 10 series of produce temperature measurements at different locations throughout the bulk. The graph shows an instable temperature simulation due to an underestimated fan capacity (m3/h) during external ventilation, in combination with the modelled temperature dependence of respiration. Therefore, model fitting was performed using the search routine lsqnonlin in Matlab. The tuning parameters were selected on the basis of physical insight, knowledge of the existing sources of model uncertainty, and the observed prediction errors. When tuning the parameters p1,resp, p2,resp and p3,resp, very small values for p2,resp and

p3,resp were found. Therefore the model was modified by using p2,resp ¼ p3,resp ¼ 0 and fitting p1,resp to the data, resulting in an improved model fit. Fitting the fan capacity (m3/h) as well yields a tuned fan capacity of a factor 1.9 larger than the specified (theoretical) fan capacity. Also, heat transfer from outside the lock to inside the lock was influenced more by the ambient temperature and less by the ground temperature than anticipated on theoretical grounds. The constant ground temperature determined only 45% of the heat transfer through the floor in the lock. Using the five new parameters that take these effects into account, the model predicted the mean produce temperature perfectly well (top right graph of Fig. 4). The bottom part of Fig. 4 depicts cross validations of the model on data of a different store in the same season and of the same store in the previous season. Both showed good results without any additional change in parameters and hence the model’s ability to predict mean produce temperatures is validated. Figs. 5 and 6 depict some produce temperatures predicted by a 10-layer model with fitted parameters. Fig. 5 depicts the merely physics-based prediction of produce temperature pull-down curves at different heights in the bulk, when ventilated with ambient air of 6  C, 0% RH. It is interesting to observe that temperature in the bulk drops below 6  C due validation case 1

case 1: original parameters 16

measurement

14

20

model

Tpr (ºC)

Tpr (ºC)

12 15

10 8

10

6 5 0

50

100

150

4

200

0

50

time (days)

100

150

200

time (days)

validation case 2

validation case 1, previous season

12

10 9

10

Tpr (ºC)

Tpr (ºC)

8 8

7 6

6 5 4

0

50

100

time (days)

150

200

4

0

10

20

30

40

time (days)

Fig. 4. Modelled and measured bulk temperatures for calibration and validation of the model.

50

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204

12

step response 10−layer and 1−layer model, input air: T = 6 ºC, RH = 0% Text Tpr 1−layer

11

Tpr, 1−10 10−layer

Tpr (ºC)

10 9

203

The resulting calibrated model accurately predicted the measured produce temperature throughout a complete storage season for three separate cases (Fig. 4), while only five calibration parameters had been calibrated by fitting the model to one of the three datasets. The model’s sound physical basis, the good data-fit and the orderly documentation assure the suitability of the resulting bulk-store model as a tool to quickly predict the consequences of specific choices in design and/or control parameters.

8

Acknowledgements 7

to evaporative cooling. This is a known phenomenon in practice. In Fig. 6 the top and bottom layer produce temperatures are shown together with the produce temperature predicted by the one-layer model.

This work was part of the WIC2 project (www. afsg.nl/wic2), which was carried out by the co-operating parties Wageningen UR, PRIVA, Tolsma Techniek Emmeloord and Weathernews International. We are especially indebted to Tolsma Techniek Emmeloord for sharing the measurement data with us. This project was supported with a grant of the Dutch Programme EET (Economy, Ecology, Technology) a joint initiative of the Ministries of Economic Affairs, Education, Culture and Sciences and of Housing, Spatial Planning and the Environment. The programme is run by the EET Programme Office, SenterNovem.

4. Conclusions

References

6 5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

time (days) Fig. 5. Temperature step response for one and ten layer model.

A modelling procedure was presented which combines the zonal decomposition presented by [1] with a processoriented matrixevector notation [2]. The procedure’s suitability to model distributed climate dynamics in storage facilities was illustrated by its successful application to the case of potato storage. The resulting model presentation facilitates an easy understanding of the mutual interactions between processes and states.

16

case 1: calibrated parameters, simulation results for one and 10 layer models model one layer model 10 layers: bottom model 10 layers: top

14

Tpr (ºC)

12

10

8

6

4

0

20

40

60

80

100

120

140

160

180

200

time (days) Fig. 6. Predicted bulk temperatures for one and ten layer model.

[1] D.J. Tanner, A.C. Cleland, L.U. Opara, T.R. Robertson, A generalised mathematical modelling methodology for design of horticultural food packages exposed to refrigerated conditions: Part 1, formulation, International Journal of Refrigeration 25 (2002) 33e42. [2] M. Henze, C.P.L. Grady Jr., W. Gujer, G.V.R. Marais, T. Matsuo, Activated Sludge Model No. 1. IAWQ Scientific and Technical Report no. 1, IAWQ, London, UK, 1987. [3] F.W. Bakker-Arkema, Selected aspects of crop processing and storage: a review, Journal of Agricultural Engineering Research 30 (1984) 1e22. [4] K. Gottschalk, Mathematical modelling of the thermal behaviour of stored potatoes & developing of fuzzy control algorithms to optimise the climate in storehouses, in: W. Day, P.C. Young (Eds.), Proc. Math. and Control Appl. in Agric. and Hortic., Acta Horticultura 406, 1996, pp. 331e339. [5] K.K. Khankari, Application of a numerical model for prediction of moisture migration in stored grain, Transactions of the ASAE 38 (6) (1995) 1789e1804. [6] A.N. Marchant, P.H. Lidstone, T.W. Davies, Artificial intelligence techniques for the control of refrigerated potato stores. Part 2: heat and mass transfer simulation, Journal of Agricultural Engineering Research 58 (1994) 27e36. [7] L.A. Schaper, Potato storage ventilation duct simulation on a PC, Applied Engineering in Agriculture 7 (4) (1991) 465e472. [8] Y. Xu, D. Burfoot, Predicting condensation in bulks of foodstuffs, Journal of Food Engineering 40 (1999) 121e127. [9] H.B. Nahor, N. Scheerlinck, P. Verboven, J. van Impe, B.M. Nicola€ı, Continuous/discrete simulation of controlled atmosphere (CA) cool storage systems: model development and

204

[10]

[11]

[12] [13]

L.J.S. Lukasse et al. / International Journal of Refrigeration 30 (2007) 195e204 validation using pilot plant CA cool storage, International Journal of Refrigeration 27 (2004) 884e894. H.B. Nahor, N. Scheerlinck, P. Verboven, J. van Impe, B.M. Nicola€ı, A continuous/discrete simulation of controlled atmosphere (CA) cool storage systems: validation using industrial CA cool storage, International Journal of Refrigeration 28 (2005) 461e470. A. Rastovski, A. van Es, Storage of Potatoes e Post-Harvest Behaviour, Store Design, Storage Practice, Handling (Pudoc, Wageningen, The Netherlands) (1987). C.F.H. Bishop, W.F. Maunder, Potato Mechanisation and Storage, Farming Press, 1980, 208 pp. R.G. van der Sman, Simple model for estimating heat and mass transfer in regular-shaped foods, Journal of Food Engineering 60 (4) (2003) 383e390.

[14] J.L. Woods, Moisture loss from fruits and vegetables, Postharvest-News Information 1 (3) (1990) 195e199. [15] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, Wiley, New York, 1960. [16] L.E. Gra¨hs, B. Hylmo¨, A. Johansson, C. Wikberg, The two point temperature measurement e a method to determine the rate of respiration in a potato pile, Acta Agriculturae Scandinavica 28 (1978) 231e236. [17] P.A. Schippers, The rate of respiration of potato tubers during storage. 2. Results of experiments in 1972 and 1973, Potato Research 20 (1977) 189e206. [18] W.G. Burton, Physiological and biochemical changes in the tuber as affected by storage conditions, in: Proceedings of the Fifth Triennial Conference Eur. Ass. Potato Res., Norwich, England, 1972, pp. 63e81.