Modelling the Transient Thermal Behaviour of Sand Substrate heated by Electric Cables

ARTICLE IN PRESS Biosystems Engineering (2005) 90 (2), 203–215 doi:10.1016/j.biosystemseng.2004.11.005 SE—Structures and Environment

Modelling the Transient Thermal Behaviour of Sand Substrate heated by Electric Cables M.D. Ferna´ndez1; M.R. Rodrı´ guez1; F. Maseda1; R. Velo1; M.A. Gonza´lez1 1

Department of Agroforestry Engineering, Escola Polite´cnica Superior, University of Santiago de Compostela, Campus Universitario s/n, 27002 Lugo, Spain; e-mail of corresponding author: [email protected] (Received 16 February 2004; accepted in revised form 18 November 2004; published online 26 January 2005)

The development of optimal control strategies in greenhouses requires precise models of energy transfer. A flexible model that enables the introduction of different substrate parameters and boundary conditions was developed for electric cable heating installations. The model is based on finite element analysis and is implemented by using a commercial code, which provides the two-dimensional description of the thermal state of heated greenhouse substrates. The model presented includes methods for the estimation of apparent thermal conductivity and methods for the estimation of heat capacity per volume given other physical properties of the substrate. The model allows for the variations of moisture and compaction in different substrate layers. The final model is compared with the experimental data obtained from an electric cable heating installation applied to a bench. Temperature data was measured at six locations of the substrate volume. A root mean square error (RMSE) of 0
56 1C was obtained after implementation during a 3-day simulation. This model can be used to monitor and evaluate substrate heating systems. r 2004 Silsoe Research Institute. All rights reserved Published by Elsevier Ltd

1. Introduction Greenhouse conditioning is oriented to the application of optimal control strategies. It is mainly based on obtaining the highest performance with minimum resource consumption. Special attention is given to the energy efficiency of installations. This global concept has an important application in heating systems, and includes automation of the installation and monitoring of its effects. Heating systems represent the highest consumption of energy in greenhouses (Hanan, 1998). Precise models of the behaviour of the greenhouse components that have an influence on plant development and growth are required to implement optimal control. The present paper focuses on electric cable heating in substrates, which favours crop growth and development (Washitani & Saeki, 1986; Andales et al., 2000; Wurr et al., 2001). This system shows good properties for efficient regulation and control, although it demands higher energy costs. Yet, the correct modelling of the heat transfer process applied to the electric cable helps to obtain optimal conditions. 1537-5110/$30.00

Many existing models are designed to predict soil temperatures under natural conditions, among these Novak (1993) and Moroizumi and Horino (2002). These models do not enable the introduction of the heat supplied by the heating system, and they do not consider the effects of the heat supplied. Therefore, a model that can be applied under conditions of artificial substrate heating is needed, in order to improve the environmental conditions for crop development. Such a model can be distinguished from models subject to natural conditions, in four ways. Firstly, the geometry of the system to be modelled is more complex due to the heat supplied by the heating system and its distribution within the substrate. Secondly, the thermal properties of heated substrates show greater variability than the thermal properties of substrates under environmental conditions. These depend on other physical properties such as texture, degree of compaction or water content (de Vries, 1963; Ochsner et al., 2001). The conditions that affect the thermal properties can be dependent on two variables: space and time. Also, studies by other authors (Campbell et al., 203

r 2004 Silsoe Research Institute. All rights reserved Published by Elsevier Ltd

ARTICLE IN PRESS 204

M.D. FERNA´NDEZ ET AL.

Notation

A a1 a2 [B] Cp Cv D Db Dbd Ds ERMS e h hf

K [K] Ke [Ke] Kxx, Kyy fL g L1

L2 L3

surface of the element, m2 depth of the lowest substrate layer, mm depth of the second substrate layer, mm matrix resulting from ½B ¼ fLgfNgT heat capacity per mass unit, J kg1 1C1 volumetric heat capacity, J m3 1C1 density, kg m3 bulk density, kg m3 bulk density of dry substrate, kg m3 solid density, kg m3 root mean square error, 1C insulator thickness, mm height of the insulator over the substrate, mm surface convection coefficient evaluated at a mean temperature between both sides of the surface, W m2 thermal conductivity, W m1 1C1 thermal conductivity matrix substrate effective thermal conductivity, W m1 1C1 effective thermal conductivity matrix thermal conductivity of the element in the abscissa and ordinate directions, respectively, W m1 1C1 Laplace transform line in a two-dimensional space, representative of a surface on which a certain temperature acts line in a two-dimensional space, representative of a surface on which a heat flow acts line in a two-dimensional space, representative of a surface on which a temperature acts as surface convection

1994; Tarnawski et al., 2000) report a variation in the apparent thermal conductivity depending on temperature. Considering the thermal gradients produced within heated substrates, the variation of apparent thermal conductivity cannot be disregarded in thermal analysis. Thirdly, the temporal variation of temperatures requires the thermal analysis to be performed in a transient state, which allows modelling for long periods of time, determined by the crop growth cycle. It also allows modelling at short-time intervals, which can be applied as a base for control (Challa & van Straten, 1993). Finally, the generated model has to be accessible to potential users, mainly, greenhouse crop producers, technicians, and suppliers of heating installations, by

N fng fqg q rH s0 s1 T T TB Tc fT e g Tm TS t t1 t2 t3 fvg dT y fZg r r f gT

number of observations set of element shape functions heat flow vector, W m2 heat flow acting on a surface, W m2 rate of generation of heat per volume unit, W m3 additional spacing to minimum initial spacing s1, mm semi-spacing between cables, mm temperature, 1C temperature acting on a surface, 1C temperature at the surface of convection, 1C calculated temperature, 1C vector of the nodal temperature of an element, 1C measured temperature, 1C temperature in the fluid adjacent to the surface, 1C time, s depth of the third substrate layer, mm depth of the fourth substrate layer, mm depth of the surface substrate layer, mm mass transport velocity vector, m s1 virtual variation of temperature, 1C substrate moisture measured over the total volume, m3 m3 normal unitary vector to the surface grad operator divergence operator transpose of a matrix

using a commercial code (ANSYS v.6.0). The generated model enables the user to perform simulations by modifying parameters related to substrate, geometry of the heating system, and environmental conditions, without affecting the other features of programming. To model the physical phenomena occurring in greenhouses, resolutions based on analytical methods were used (A´lvarez et al., 1996). However, there is a tendency towards the use of numerical methods (Tavares et al., 2001), such as finite differences (de la Plaza et al., 1999); electrical analogy (Guaraglia et al., 2001); neural networks (George, 2001; Mihalakokou, 2002); and finite elements (Puri, 1987). This tendency responds to the limitations of analytical methods as regards modelling requirements. In this case, a

ARTICLE IN PRESS 205

MODELLING THE TRANSIENT THERMAL BEHAVIOUR

commercial code for finite element analysis was used. Yet, all the presented models contain a series of substrate simplifications that restrict their precision and descriptive level. Some of the models presented are one-dimensional (A´lvarez et al., 1996; Guaraglia et al., 2001; de la Plaza et al., 1999) and therefore, show limitations where heat transfer is two-dimensional. Although a two-dimensional model was proposed by Kurpaska and Slipek (2000), its purpose was not to represent the evolution of temperatures occurring within the studied section, but simply to optimise the heating system. For this reason, modelling was performed taking into consideration a steady state. Other simplifications present in the existing models refer to the characterisation of the substrate. The purpose of this paper is to develop a twodimensional model, based on finite elements, capable of describing the thermal state of substrates heated by means of electric cable heating, given the geometry of the heating system and the thermal properties of the substrate and the insulator. Unlike the previous models (A´lvarez et al., 1996; Kurpaska and Slipek, 2000; de la Plaza et al., 1999), this model enables: (1) the introduction of different properties of compaction and moisture, and mixture of materials in the modelled substrate; (2) the transient analysis of the system by showing at each moment the two-dimensional distribution of temperatures, and by allowing the assessment of heat flows; and (3) the use of a minimum number of boundary conditions—initial temperature of substrate and temperatures at the surface of the substrate and at the heating cable at each step of analysis. Considered in this manner, the model can be incorporated into optimal control strategies, permitting: (1) the monitoring of the thermal state of heated greenhouse substrates with a minimum number of known variables; and (2) the assessment of substrate heating systems regarding energy and heating quality. Such assessment enables comparison of different heating systems under suitable environmental and control conditions without threatening greenhouse production.

2. Theoretical considerations In this study, heat flow is considered two-dimensional and heat transfer is restricted to conduction and convection. Radiation can be neglected (Buchan, 1991; Jury et al., 1991). The expressions defining heat flow are ruled by the first law of thermodynamics of the conservation of energy. When referred to a differential control volume, represented in the two-dimensional

model as an area, it is written as
 qT DC p þ fvgT fLgT þ ðfLgT fqgÞ ¼ rH qt

(1)

where: D is density in kg m3; Cp is heat capacity per mass unit in J kg1 1C1; T is temperature in 1C; and t is time in s; fLg is Laplace transform, in a twodimensional space and expressed as a matrix as (q) fL g ¼

qx q qy

(2)

and fvg is the vector of mass transport velocity in m s1, developed within a two-dimensional space as ( ) vx fvg ¼ (3) vy fqg is the heat flow vector in W m2; rH is the rate of generation of heat per volume unit in W m3; and the superscript T means the transpose of a matrix. The terms fLgT and fLgT fqg may also be interpreted as rT and r
q; where r represents the grad operator and r
the divergence operator. The mechanisms of heat transfer in the insulator are conduction and convection, while in the substrate the only mechanism taken into consideration is transfer by conduction. Fourier’s law was used in several of the physical models proposed (Persaud & Chang, 1983; Novak, 1986; Karam, 2000) to relate the heat flow vector to the thermal gradients in the study of heat transfer in soils. In order to use this law within the finite element method, it must be expressed as a matrix fqg ¼ ½K fLgT

(4)

where [K] is the thermal conductivity matrix, which can be developed within a two-dimensional space as " # K xx 0 ½K  ¼ (5) 0 K yy where Kxx and Kyy stand for the thermal conductivity of the element in the abscissa and ordinate directions, respectively. The effective thermal conductivity matrix [Ke], in W m1 1C1, is used for soil instead of [K], and includes the effects of conduction and convection in the liquid and gaseous phases of soil (Jury et al., 1991). Fourier’s law can therefore be used as a description of coupled flow of heat and moisture, in which the flow conducted by water vapour is considered, and the heat transferred by liquid water is neglected (Buchan, 1991). Combining Eqn (1) and Eqn (4), the following expression is obtained:
 qT DC p þ fvgT fLgT ¼ fLgT ð½K fLgT Þ þ rH (6) qt

ARTICLE IN PRESS M.D. FERNA´NDEZ ET AL.

206

Expanding Eqn (6) within a two-dimensional coordinate system, the following expression is obtained:
 qT qT qT þ vx þ vy DC p qt qx qy

 q qT q qT K xx K yy ¼ rH þ þ ð7Þ qx qx qy qy The code used, ANSYS v.6.0, enables the application of three types of boundary conditions. It is assumed that these conditions cover the entire element. (1) Specification of the temperatures acting on a surface is represented by a line L1 in a two-dimensional space: T ¼T

n

(8)

where T* is the specified temperature in 1C. (2) Specification of the heat flows acting on a surface is represented by the linear element L2 in the twodimensional model: fqgT fgg ¼ q

(9)

where: fgg is the normal unitary vector to the surface (in a three-dimensional coordinate system) or to the line (in a two-dimensional coordinate system); and qn is the specified heat flow, in W m2. (3) Specification of convection on surfaces is represented in a two-dimensional space by lines L3. These are ruled by Newton’s law of cooling: fqgT fgg ¼ hf ðT B  T S Þ

(10)

where: hf is the surface convection coefficient, evaluated at an average temperature between both sides of the surface, in W m2; TB represents the temperature in the fluid adjacent to the surface, in 1C, and TS is the temperature in the surface of convection, in 1C. The derivation of heat flow matrices applied by the finite element code ANSYS v.6.0 (ANSYS, 2001) is presented below. Combining Eqns (4), (9) and (10), the following equations are obtained: fggT ½K fLgT ¼ qn

(11)

fggT ½K fLgT ¼ hf ðT B  TÞ

(12)

Multiplying Eqn (6) by a virtual variation of temperature dT in 1C, which can vary according to position and time of analysis; integrating the Eqn (6) over the volume of an element represented by a surface A in m2, in a two-dimensional coordinate system; and combining Eqns (11) and (12), the following expression

is obtained:
 Z
qT þ fvgT fLgT DC p dT qt A  T þ fLg ðdT Þ½K fLgT dðAÞ Z Z ¼ dTqn dðL2 Þ þ dThf ðT B  TÞ dðL3 Þ L2 L3 Z dTrH dðAÞ þ

ð13Þ

A

In order to develop the matrix equations of each element, T can vary both in time and in space, and so can dT: This dependence is established as T ¼ fngT fT e g

(14)

where: fng ¼ fnðx; y; tÞg is the set of element shape functions, also as a function of space and time, and fT e g ¼ fT e ðTÞg is the vector of the nodal temperature of an element, in 1C. Deriving Eqn (14) in terms of time, the following expression is obtained:

qT T qT e ¼ fng (15) qt qt The virtual variation of temperature dT is expressed in the same way dT ¼ fdT e gT fng

(16)

The term fLgT is written as fLgT ¼ ½BfT e g

(17)

where [B] is the matrix resulting from ½B ¼ fLgfngT

(18)

Thus, the variational state of Eqn (13) can be combined with Eqns (14)–(18), in order to obtain

Z T T qT e DC p fdT e g fngfng dðAÞ qt A Z DC p fdT e gT fngfvgT ½BfT e g dðAÞ þ ZA þ fdT e gT ½B T ½K ½B fT e g dðAÞ A Z fdT e gT fngqn dðL2 Þ ¼ L2 Z þ fdT e gT fnghf ðT B  fngT fT e gÞ dðL3 Þ ZL3 ð19Þ þ fdT e gT fngrH dðAÞ A

In Eqn (19), the density of the material is assumed to remain constant over the volume of the element. On the other hand, the specific heat and the rate of generation of heat per volume unit can vary within each element. Finally, fT e g; fqT e =qtg and fdT e g are nodal quantities

ARTICLE IN PRESS 207

MODELLING THE TRANSIENT THERMAL BEHAVIOUR

Sandwich panel compartment

Position T1

Heating cable

50

50

7575 75 75

300

Y

T2

225

X

75 150 300

275

T1 00 75 150 225 300

250

H 75 150 225

250 Section D - D'

Position T1 50

Position T2

Position H 50

1200

3. Materials and methods

D'

D

The step of analysis must be selected to allow for the convergence of the solution. The boundary conditions introduced at each step of analysis are representative of the variations experienced by real continuous variables. Equation (20) defines the thermal behaviour of each of the elements that form the mesh that separates substrate from insulator. The application of the equation requires the introduction of the following properties for the materials: D, K and Cp. These properties are constant and known in the case of the insulator, but not in the case of the substrate. In the case of the substrate, the following thermal properties must be calculated: (1) bulk density Db, used instead of D; (2) Ke, used instead of K, and (3) volumetric heat capacity Cv, related to Cp by bulk density. These properties are calculated from texture data, solid density Ds, bulk density of dry substrate Dbd and substrate moisture y:

1600

A

1700

L3

Position H Position T2

600

which do not vary on the element, and so can also be eliminated from the integral. Because all the quantities were previously multiplied by a vector fdT e g; this term can be dropped from the equation

Z qT e D C p fngfngT dðAÞ qt A Z þ D C p fngfvgT ½B dðAÞfT e g Z A ½B T ½K ½B  dðAÞfT e g þ A Z Z ¼ fngqn dðL2 Þ þ T B hf fng dðL3 Þ L2 L3 Z Z T hf fngfng fT e g dðL3 Þ þ rH fng dðAÞ ð20Þ 

Plan

Fig. 1. Components and arrangement of the substrate heating cable installation. Position of temperature probes: T1, at the vertical of the electric cable; T2, between cables and H, position of humidity probes

3.1. Experimental test An experimental test was required to develop a model allowing the thermal analysis of substrate heated by means of electric cable heating. The devised test made it possible to obtain the boundary conditions to apply to the model, as well as data of substrate moisture at different depths and, additionally, temperature data to validate the model. The experimental test was conducted in the experimental greenhouse of the Department of Plant Production of the University of Santiago de Compostela, located in Lugo, Spain. An electric substrate heating installation was implemented. The heating installation consisted of a compartment built with a sandwich panel structure with a core of extruded polystyrene foam. Figure 1 shows the dimensions of the compartment.

Inside, a 230 V AC-fed parallel heating cable, 30 W m1, with silicone insulation (Model AKO-5234, AKO, Barcelona, Spain) was arranged and actuated by means of a contactor. The volumetric water content of substrate was maintained at a range of 0
10–0
14 by means of an irrigation system. The irrigation system was made up of 12 microsprinklers, 20 l h1, with an effective diameter of 1200 mm, working pressure of 200 kPa, located at a height of 300 mm over the substrate surface. It was opened and closed through a solenoid valve, 28 m3 h1, 24 V DC, with an operating pressure ranging from 100–600 kPa, which agreed with data obtained from the humidity probes. The contacts of the solenoid valve opened when subjected to pressure. Inside each compartment, sand substrate with a 379 mm mean geometric diameter was arranged. The

ARTICLE IN PRESS M.D. FERNA´NDEZ ET AL.

e

s1

s0

t3

Possible location of heating cable

h

solid density of substrate was 2686 kg m3, and the average bulk density of dry substrate was 1256 kg m3, variable with depth. Negative thermal coefficient thermistor temperature sensors (Models Campbell 107 & Campbell 108, Campbell Sci. Ltd, Loughborough, UK) and matric potential sensors (Model Watermark 257, Delta-T Devices Ltd, Cambridge, UK), with porous matrix and measurement range 0–200 kPa and 5 V AC excitation were placed at different substrate locations. Sensors enabled data related to these variables to be collected for the operation of the heating and irrigation installations. Figure 1 shows the location of sensors and their designation. Measurements were taken by temperature and humidity probes at 1-min intervals. The mean values were stored at 15 min intervals by means of a data logger (Model CR10X, Campbell Sci. Ltd, Loughborough, UK). The data logger was connected to a programmable controller (Model S7-200, Siemens AG, Mu¨nchen, Germany), which controlled the installations given the received data.

e a1 a2 t1 t2

208

(a)

(b)

3.2. Formulation of the model according to the finite element method The need to develop an accessible model based on the finite element method, determined the use of commercial codes. The modelling conducted was framed within complex geometry transient thermal analysis, with a mixture of materials and a change of physical thermal properties depending on temperature and time. In order to meet the mentioned requirements, the commercial code ANSYS v.6.0 (Ansys Inc., Canonsburg, PA) was selected for use. The accessibility of the code was determinant for its selection. 3.2.1. Pre-process The development of the model was structured in several phases. During the first (pre-processing) phase the geometrical definition and the discretisation of the problem were performed. The model must be capable of simulating all the conditions of operation of a substrate heating system. Therefore, the model must be able to accommodate variations in bench dimensions and heating cable positions, both in depth and in spacing. For simplicity, these variables were parameterised in the model as shown in Fig. 2(a). Furthermore, the two-dimensional model permitted the introduction of the physical properties of the materials at different depths. For that purpose, substrate was divided into layers. Once the geometry was defined, the geometrical system proposed was meshed to obtain the elements to which the function defined in Eqn (20) was applied. The

Fig. 2. (a) Parameterised geometry developed for the finite element model. Parameters: e, insulator thickness; s0, additional spacing to minimum initial spacing s1; s1, semi-spacing between cables; a1, depth of the lowest substrate layer, a2, depth of the second substrate layer; t1, depth of the third substrate layer; t2, depth of the fourth substrate layer; t3, depth of the surface substrate layer; h, height of the insulator over the substrate, and (b) mesh

mesh used to discretise the model was rectangular and fine [Fig. 2(b)]. The average size of the elements was 100 mm2. The elements were smaller in the zones where the boundary conditions were applied (location of the heating cable and surface of the substrate) and larger in intermediate areas. A shape function for the physical behaviour of each of the elements constituting the model was included. The shape function was defined by Eqn (20), and it was also defined by the ANSYS v.6.0 element PLANE55. The global stiffness matrix for the whole set of elements in the model was built by reassembling the set of equations. 3.2.2. Solution In the second phase, solution, the stiffness matrix was solved simultaneously for all the elements of the mesh, defined during the pre-process. The results were nodal values (temperature) and elemental values (heat flow). The solution method used was sparse direct solver, which provided a good rate of convergence and robustness of analysis. During the solution phase, the thermal loads were applied and results were obtained. Thermal loads were fixed parametrically without affecting the rest of the

ARTICLE IN PRESS MODELLING THE TRANSIENT THERMAL BEHAVIOUR

analysis. The introduction of the initial boundary conditions in a thermal analysis required the establishment of all the nodal temperatures. If the boundary conditions were initially unknown, a uniform initial temperature was assumed for all the nodes. As a result, an initial period occurred during which results were affected by an ideal assumption of the initial conditions. Each analysis was made up of 97 steps (a first step to determine the initial conditions and 96 further steps) that allowed the study of 1 day, divided into steps of 15 min each. The following boundary conditions were applied to each step of analysis, made up of 1 min long substeps: temperature at the surface of the substrate, air temperature at 1 m over the surface of the substrate, expressed as the surface convection at the surface of the insulator, and temperature of the electric cable. The loads were considered as ramped loads, between the steps of analysis. Newton–Raphson algorithm was used to obtain the results. The thermal properties were evaluated at each substep of analysis, with a maximum of 15 iterations. The solution was achieved when the difference between the results obtained in two subsequent iterations did not exceed tolerance values. The following tolerance values were used: 0
001 1C for temperatures, and 106 W m2 for the heat balance of the whole system. 3.2.3. Post-process The large quantity of results produced by the numerical method demanded the completion of a last phase, the post-process, in order to obtain the required results. In this case, results were required to conduct the experimental validation of the model, to further develop the analysis during the following days, and to describe the physical phenomenon modelled. 3.3. Validation The results obtained by the final model were compared with experimental data by applying root mean square error (RMSE): !0:5 N 1X 2 E RMS ¼ ðT c  T m Þ (21) N 1 where: ERMS is the root mean square error; N is the number of observations; Tc is the calculated temperature; and Tm is the measured temperature.

4. Results and discussion 4.1. Development of the models After the formulation of the model based on finite elements, the model must be developed. Development

209

involves determining the degree of complexity of the model as regards the physical properties of soil and their variation. Table 1 shows the considerations taken to develop the model until the last step, which constitutes the final model used. In order to analyse the effects of the environmental conditions and the effects of artificial heat supply, a day is selected to start running the heating system. During the initial phase, the temperature of the whole substrate is considered uniform and estimated as an average of the temperatures measured at the initial stage of the experimental test (Fig. 3). The boundary conditions, applied at 15 min intervals, correspond to the measurements of the experimental test conducted on 28th January 2002 (Fig. 4). The results obtained for each assumption are composed of nodal temperatures and heat flows at each substep of analysis. The comparison between the simulated situation and the experimental situation is conducted by means of temporal graphs for temperature at positions T1 and T2 (Fig. 1) in models A-1 (Fig. 4) and A-4 (Fig. 5). In model A-1 (Fig. 4), the measured and the simulated temperatures overlap until the artificial heating source constituted by the heating system starts running (at a time of 48 600 s). At this point, the differences between the measured and the simulated temperatures become evident. The simulated temperatures are shown to be lower than the measured temperatures. With increasing temperature, the difference between both temperatures becomes more pronounced, reaching a value of 8
83 1C (Table 2). In addition, the increase in simulated temperature is slower than the increase in temperature measured in the experimental test; a delay in the response of the model occurs. Taking into consideration that the tendency shown by the model is adequate, it is considered necessary to evaluate the variation of effective thermal conductivity Ke as a function of temperature, instead of using a constant conductivity obtained experimentally (Table 1). That variation is introduced in model A-2 (Table 1) by using the method by de Vries (1963), modified by de la Plaza et al. (1999). As a result, model A-2 shows differences between measured and simulated temperatures that do not reach values of 3 1C (Table 2), clearly improving previous results. Although this model provides a quicker response of the simulated temperatures than the previous model, delay still occurs. The results of this model (A-2) show how the introduction of the variation of effective thermal conductivity depending on temperature constitutes a crucial factor in the analysis of heat transfer in heated substrates. However, the temperatures recorded at depth are too low, and the temperatures recorded at the surface are too high. This

ARTICLE IN PRESS M.D. FERNA´NDEZ ET AL.

210

Table 1 Considerations taken throughout the development of the model and properties of the substrate under a solid density of 2686 kg m3, and the following texture: 96.15% sand, 0.11% silt, 3.74% clay, and 1.13% organic matter Model

Depth, mm

Physical properties Bulk density of dry substrate (Dbd), kg m3

Determination of substrate thermal properties

Water content (y), m3 m3

Volumetric heat capacity (Cv), kJ m3

Effective thermal conductivity (Ke), W m1 1C1 at temperatures of 10 1C

A-1

0–300

A-2

0–75 75–150 150–225 225–300

1256 1256

30 1C

60 1C

0
368a

0
133

922
63

0
139 0
134 0
129 0
129

955
43 936
05 916
14 915
61

0
987b 0
972b 0
956b 0
955b

1
029b 1
013b 0
999b 0
997b

1
216b 1
200b 1
190b 1
190b

A-3

0–75 75–150 150–225 225–300

1172 1230 1300 1360

0
139 0
134 0
129 0
129

955
43 936
05 916
14 915
61

0
896b 0
943b 1
005b 1
078b

0
935b 0
984b 1
048b 1
123b

1
112b 1
171b 1
247b 1
331b

A-4

0–75 75–150 150–225 225–300

1172 1230 1300 1360

0
139 0
134 0
129 0
129

955
43 936
05 916
14 915
61

0
982c 0
981c 1
046c 1
121c

1
026c 1
025c 1
092c 1
170c

1
222c 1
226c 1
307c 1
396c

A-5

0–10 10–75 75–150 150–225 225–300

1172 1172 1230 1300 1360

0
070 0
139 0
134 0
129 0
129

772
73 955
43 936
05 916
14 915
61

0
737c 0
982c 0
981c 1
046c 1
121c

0
774c 1
026c 1
025c 1
092c 1
170c

1
018c 1
222c 1
226c 1
307c 1
396c

a

Experimental value obtained by using the method of Kaminsky (1994). Values calculated by applying the method of de Vries (1963) modified according to de la Plaza et al. (1999). c Values calculated by applying the method of Campbell et al. (1994). b

Tamb

T1-00-M

T1-225-M

Temperature, ºC

35

25

15

5 0

14 400

28 800

43 200 Time, s

57 600

72 000

86 400

Fig. 3. Evolution of temperatures used as boundary conditions: ambient-air temperature in the greenhouse, at the surface of substrate and at the electric cable (225 mm depth) during the experimental test on 28th January 2002

contrast shows the heterogeneity of the physical properties of soil based on depth. Model A-3 incorporates the variation of bulk density Db (Table 1) in four layers (75 mm thick) of the

substrate. A slight improvement is achieved for all the measured temperatures, as compared with the assumption of constant bulk density (Table 2). Considering the results offered by the previous models, it can be inferred that the major problem posed by thermal properties in heated substrates and, particularly, by effective thermal conductivity, is its variation depending on temperature. This variation was modelled, among other authors, by Campbell et al. (1994). This model enables the estimation of the thermal conductivity of the liquid and gaseous phases of a soil or substrate at different temperatures as a continuous function. Few empirical parameters are required to formulate this model. The results for model A-4 (Fig. 5) show a decrease in the differences between measured and simulated temperatures. This decrease makes the method by Campbell et al. (1994) more appropriate for temperature simulation in heated substrates. Many authors have verified the great influence of the surface layer on soil heat flow (Sikora et al., 1990; Luo et al., 1992; Krarti et al., 1995). The characteristics of the surface layer differ from the characteristics of the rest of the substrate. One of these characteristics is soil

ARTICLE IN PRESS 211

MODELLING THE TRANSIENT THERMAL BEHAVIOUR

T1-75-M

T1-150-M

T1-300-M

T1-75-S

T1-150-S

T1-300-S

T1-150-M

T1-300-M

T1-75-S

T1-150-S

T1-300-S

35

25

Temperature, ºC

Temperature, ºC

35

T1-75-M

15

25

15

5 0

14 400

28 800

(a)

43 200 57 600 Time, s

72 000

5

86 400

0

14 400

28 800

(a)

T2-75-M

T2-150-M

T2-300-M

T2-75-S

T2-150-S

T2-300-S

35

43 200 Time, s

57 600

72 000

T2-75-M

T2-150-M

T2-300-M

T2-75-S

T2-150-S

T2-300-S

86 400

25 Temperature, ºC

Temperature, ºC

35

15

25

15

5 0 (b)

14 400

28 800

43 200 Time, s

57 600

72 000

86 400

Fig. 4. Evolution of temperatures at depths of 75, 150, and 300 mm, measured (-M), and simulated (-S) by model A-1, at positions (a) T1 and (b) T2, during the experimental test on 28th January 2002

moisture content. Normally, the surface layer presents a lower moisture content than the rest of the substrate, except with rain fall or irrigation. Therefore, it shows lower thermal conductivity and lower specific heat. Model A-5 enables a more detailed analysis of the surface layer (Table 1). In the case of the sand substrate used, the thickness of the surface layer is fixed at 10 mm, after measurement during the experimental test. The humidity assigned to this layer is obtained as a linear approximation between the moisture value of the underlying layer and a zero value corresponding to the moisture content at the substrate surface, which is in contact with the environment. The results obtained by model A-5 show a slight improvement as compared with model A-4 (Table 2). In this analysis, the presence of this layer particularly affects those situations under heat supply from environ-

5 0

(b)

14 400

28 800

43 200

57 600

72 000

86 400

Time, s

Fig. 5. Evolution of temperatures at depths of 75, 150, and 300 mm, measured (-M), and simulated (-S) by model A-4, at positions (a) T1 and (b) T2, during the experimental test on 28th January 2002

ment to substrate. Under these circumstances, model A4 shows the lowest precision. Model A-5 is preferred because of the special characteristics that the surface layer shows in substrates. In spite of the increasing adjustment to reality observed throughout the evolution of the models, substantial differences persist in temperatures measured at a depth of 300 mm. The greatest differences between measured and simulated temperatures are recorded at these locations. The simulated temperatures are lower than the measured temperatures with increasing temperature, and higher with decreasing temperature. Such differences can be due to soil dessication in zones subject to high temperatures, which produces a humidity gradient at locations with the same depth. This situation introduces not only temperature-dependent variations, but also position-dependent variations on the

ARTICLE IN PRESS M.D. FERNA´NDEZ ET AL.

212

Table 2 Notation and depth of measured and simulated probes; maximum differences between measured and simulated temperatures in absolute value, and root mean square error (RMSE) of a 3-day simulation period, not taking into consideration the first 4 hours, affected by the hypothetical initial conditions Probe designation

Depth, mm

Maximum difference between measured and simulated temperature, 1C for models A-1

A-2

A-3

A-4

A-5

RMSE for A-5, 1C

T1

75 150 300

3
23 5
14 7
51

1
54 1
42 2
33

1
46 1
32 2
24

0
98 0
84 2
36

0
50 0
39 1
38

0
36 0
24 0
79

T2

75 150 300

2
70 5
07 8
83

2
58 1
56 2
32

2
38 1
44 1
76

2
08 1
13 2
57

0
76 1
07 1
81

0
34 0
46 0
87

T1, T2

75, 150, 300

8
83

2
58

2
38

2
57

1
81

0
56

thermal properties. However, this phenomenon was not considered because its effects occur at the deepest layers of the substrate, out of the volume for root growth. 4.2. Final model and validation Under environmental conditions, all the previous models calculate temperatures with the same precision. However, under artificial heat supply, model A-5 provides the best results. Therefore, model A-5 is considered the definitive model. This model, solved by using the finite element method, is described in Fig. 6. The model is implemented by applying the methods for the estimation of physical properties shown in Table 1. The following aspects were considered in the final model: (1) the calculation of effective thermal conductivity as a function of temperature was conducted by applying the method by Campbell et al. (1994); (2) heat capacity per volume was calculated by following the method proposed by de Vries (1963); (3) introduction of moisture and compaction values at each substrate layer, which slightly improves the analysed temperatures, considered as a whole [Fig. 2(a)]; (4) introduction of a differentiated 10 mm thick surface layer, with lower moisture content than deeper layers, which improves the precision of the model under high ambient-air temperatures. The physical properties of the substrate are indicated as parameters in a step previous to the formulation of the analysis according to the finite element method. This step enables the estimation of the substrate thermal properties and the introduction of the geometric

variables defining the heating system without affecting the other features of programming. The input parameters for the estimation of the substrate thermal properties are the topographic height, in m, and the physical properties of each substrate layer required by the methods by de Vries (1963) and Campbell et al. (1994): Ds, in kg m3; Db, in kg m3; y; in m3 m3; geometric mean particle diameter, in mm; and power for the recirculation function. The geometric parameters are shown in Fig. 2(a). Likewise, another module permits the introduction of the temperatures required as boundary conditions. In this module, the duration time of the steps and substeps of analysis can be modified. The model was validated by comparing the results it yielded with the experimental data collected during 3 days of operation of the heating system, between 7th February and 9th February 2002. The heating installation operated under the following programming conditions: probes at position T1 (Fig. 1), at depths of 150 and 225 mm, lower than or equal to 25 1C and lower than or equal to 55 1C, respectively. A uniform initial temperature was assumed in all the nodes. Figure 7 shows boundary conditions and Fig. 8 shows modelling during the second day, under the mentioned conditions. The comparison of the results obtained in the model with the data obtained through experimental testing during 3 days shows a root mean square error (RMSE) of 0
56 1C for the whole set of analysed probes. Table 2 shows the results for each of the six analysed probes. The results of the initial 4-h period, affected by the hypothetical initial conditions, are disregarded. The model was much more precise in the root zone, at depths below 150 mm, where the simulation reached RMSEs of 0
24 to 0
46 1C above the electric cable, and of 0
79 to 0
87 1C between cables.

ARTICLE IN PRESS 213

MODELLING THE TRANSIENT THERMAL BEHAVIOUR

Topographic height

Calculate substrate properties required for the analysis: - bulk density - volumetric heat capacity - effective thermal conductivity

Substrate physical properties: - texture - Ds, Dbd ,
- geometric mean particle diameter - powerrecirculation function

Geometrical parameters of the heating system: - e, s0, s1, a1, a2, t1, t2, t 3, h

Data input Process Results output

Create geometry

Mesh and aplication of shape function and calculated physical properties in each element Initial conditions

Ambiental temperature

Temperature at the substrate surface

Step and substep time

Application of boundary conditions to each step of analysis

Solution of stiffness matrix at each substep of analysis: - Newton-Raphson algorithm - maximum of 15 iterations - tolerance values: 0.001 ºC, 10-6Wm-2 Numerical results for validation Obtaining the required results

Application of numerical results of the last step as initial conditions for next analysis

YES

Continue analysis?

Graphs for nodal temperature and heat fluxes in the studied section at selected times

NO STOP

Fig. 6. Input, modules and output of model A-5 solved by using the finite element method, where heat capacity per volume is calculated by using the method by de Vries (1963), and effective thermal conductivity is calculated by using the method by Campbell et al. (1994). Geometrical parameters: e, insulator thickness; s0, additional spacing to minimum initial spacing s1; s1, semi-spacing between cables; a1, depth of the lowest substrate layer, a2, depth of the second substrate layer; t1, depth of the third substrate layer; t2, depth of the fourth substrate layer; t3, depth of the surface substrate layer; h, height of the insulator over the substrate. Substrate properties: Dbd, bulk density of dry substrate, kg m3; Ds, solid density, kg m3; y; substrate moisture measured over the total volume, m3 m3

The model allows the user to test different strategies for the regulation and control of the heating system (geometric parameters, operation periods, maximum

and minimum temperatures at different depths, maximum and minimum temperatures of the heating element, use of ambient-air temperature as a control

ARTICLE IN PRESS M.D. FERNA´NDEZ ET AL.

214 Tamb

T1-00-M

T1-225-M

35

T1-75-M

T1-150-M

T1-300-M

T1-75-S

T1-150-S

T1-300-S

25 Temperature, ºC

Temperature, ºC

35

15

25

15

5 0

14 400

28 800

43 200

57 600

72 000

86 400 5

Time, s

Fig. 7. Evolution of temperatures used as boundary conditions: ambient-air temperature in the greenhouse, at the surface of substrate and at the electric cable (225 mm depth) during the experimental test on 8th February 2002

0

14 400

28 800

43 200

57 600

72 000

86 400

Time, s

(a) T2-75-M

T2-150-M

T2-300-M

T2-75-S

T2-150-S

T2-300-S

parameter) by defining the boundary conditions considered appropriate. This allows for the estimation of the two-dimensional thermal state of the substrate, which enables the user to study energy efficiency, the quality of the heating system and of the control strategies applied, and to use the model to monitor heated substrates. In addition, the model enables the user to model long periods and to modify operation parameters at short intervals. The model permits modelling a crop cycle and adapting the heating system to the various stages.

5. Conclusions Given the development of the model of electric cable heating, in the case of a sand substrate used in a greenhouse, with variations of ambient-air temperature and no artificial heat supply, the finite element modelling achieves good agreement in all the points as compared to the experimental data, independently of the values and methods used for the determination of the thermal properties of substrate, and of the introduction of substrate properties based on depth. However, when an artificial heat source is present (in this case, the heating system) the precision of simulation depends on several factors. (1) The incorporation of variation in effective thermal conductivity as a function of temperature produces relevant improvements in simulation. Effective thermal conductivity suffers larger variations at high temperatures. Therefore, the use of constant values entails lower precision of simulation, mainly near the heating cable and, subsequently, in the root zone.

Temperature, ºC

35

25

15

5 0 (b)

14 400

28 800

43 200

57 600

72 000

86 400

Time, s

Fig. 8. Evolution of temperatures during the experimental test on 8th February 2002 at depths of 75, 150, and 300, measured (M) and simulated (-S) by model A-5, at positions (a) T1 and (b) T2

(2) The results of simulation vary depending on the method used to determine effective thermal conductivity. The method by Campbell et al. (1994) shows a better response than the method by de Vries (1963). (3) The consideration of variation in the properties of substrate at depth (compaction, moisture, and presence of a differentiated surface layer) affects the precision of estimation of the thermal properties of substrate and, therefore, the precision of the results of simulation. Taking into consideration the mentioned factors, a final model based on finite elements was obtained. The validation of the final model—performed by simulating three days of operation of substrate electric cable heating—yielded a root mean square error (RMSE) of 0
56 1C for the analysed set of temperatures. Higher precision was obtained above the heating cable at a depth of 150 mm, with an RMSE of 0
23 1C. The less

ARTICLE IN PRESS MODELLING THE TRANSIENT THERMAL BEHAVIOUR

precise results were obtained between cables, at 300 mm depth, with an RMSE of 0
87 1C.

References A´lvarez J; Sobro´n F; Bolado S (1996). Solution for heat flow in soil with a heat source at a fixed depth. Soil Science Society of America Journal, 60(4), 1028–1035 ANSYS (2001). ANSYS 6.0 Documentation. Ver. 6.0. ANSYS, Inc., Canonsburg Andales A A; Batchelor W D; Anderson C E (2000). Modification of a soybean model to improve soil temperature and emergence date prediction. Transactions of the ASAE, 43(1), 121–129 Buchan G D (1991). Soil temperature regime. In: Soil Analysis. Physical Methods (Smith K A; Mullins C E, eds), pp 551–612. Marcel Dekker, New York Campbell G S; Jungbauer J D; Bidlake W R; Hungerford R D (1994). Predicting the effect of temperature on soil thermalconductivity. Soil Science, 158(5), 307–313 Challa H G; van Straten (1993). Optimal diurnal climate control in greenhouses as related to greenhouse management and crop requirements. In: The Computerized Greenhouse. Automatic Control Application in Plant Production (Hashimoto Y, ed), pp 119–137. Academic Press, New York de Vries D A (1963). Thermal Properties of Soils. In: Physics of Plant Environment (van Wijk W R, ed), pp 210–235. NorthHolland, Amsterdam de la Plaza S; Benavente R; Garcı´ a J; Navas L; Luna L; Dura´n J; Retamal N (1999). Modeling an optimal design of an electric substrate heating system for greenhouse crops. Journal of Agricultural Engineering Research, 73(2), 131–139 George R K (2001). Prediction of soil temperature by using artificial neural networks algorithms. Nonlinear Analysis: Theory, Methods and Applications: An International Multidisciplinary Journal, 47(3), 1737–1748 Guaraglia D O; Pousa J L; Pilan L (2001). Predicting temperature and heat flow in a sandy soil by electrical modeling. Soil Science Society of America Journal, 65, 1074–1080 Hanan J J (1998). Greenhouses. Advanced Technology for Protected Horticulture. CRC Press, Boca Rato´n Jury W A; Gardner W R; Gardner W H (1991). Soil Physics, 5th Edn. John Wiley & Sons Inc, New York Kaminsky W (1994). Heat conduction in materials with nonhomogeneous structure. In: Experiments in Heat Transfer and Thermodynamics (Granger R A, ed), pp 19–22. Cambridge University Press, New York

215

Karam M A (2000). A thermal wave approach for heat transfer in a nonuniformal soil. Soil Science Society of America Journal, 64(4), 1219–1225 Krarti M; Claridge D; Kreider J (1995). Analytical model to predict non-homogeneous soil-temperature variation. Journal of Solar Energy Engineering—Transactions of the ASME, 117(2), 100–107 Kurpaska S; Slipek Z (2000). Optimization of greenhouse substrate heating. Journal of Agricultural Engineering Research, 76(2), 129–139 Luo Y; Loomis R T; Hsiau (1992). Simulation of soiltemperature in crops. Agricultural and Forest Meteorology, 61(1–2), 23–38 Mihalakokou G (2002). On estimating soil surface temperature profiles. Energy and Buildings, 34(3), 251–259 Moroizumi T; Horino H (2002). The effects of tillage on soil temperature and soil water. Soil Science, 167(8), 548–559 Novak M D (1986). Theoretical values of daily atmospheric and soil thermal admittances. Boundary Layer Meteorology, 34(1–2), 17–34 Novak M D (1993). Analytical solutions for two-dimensional soil heat flow with radiation surface boundary conditions. Soil Science Society of America Journal, 57(1), 30–39 Ochsner T E; Horton R; Ren T (2001). A new perspective on soil thermal properties. Soil Science Society of America Journal, 65, 1641–1647 Persaud N; Chang A C (1983). Estimating soil temperature by linear filtering of measured air temperature. Soil Science Society of America Journal, 47(5), 841–847 Puri V M (1987). Heat and mass transfer analysis and modeling in unsaturated ground soils for buried tube systems. Energy in Agriculture, 6, 179–193 Sikora E; Gupta S C; Kossowski J (1990). Soil temperature predictions from a numerical heat-flow model using variable and constant thermal diffusivities. Soil Tillage Research, 18(1), 27–36 Tarnawski V R; Leong W H; Bristow K L (2000). Developing a temperature-dependent Kersten function for soil thermal conductivity. International Journal of Energy Research, 24(15), 1335–1350 Tavares C; Gonc¸alves A; Castro P; Loureiro D; Joyce A (2001). Modeling on agriculture production greenhouses. Renewable Energy, 22(1–3), 15–20 Washitani I; Saeki T (1986). Germination responses of Pinus densiflora seeds to temperature, light and interrupted imbibition. Journal of Experimental Botany, 37, 1376–1397 Wurr D; Hanks G; Fellows J (2001). The effects of bulb storage temperature, planting date and soil temperature on the growth and development of Narcissus bulb units. The Journal of Horticultural Science and Biotechnology, 76(4), 465–473