Category Archives: EuroSun2008-7

Desiccant wheel model. Model description

The heat and mass transfer model for the desiccant wheel used below is based on the analogy method with heat transfer that occurs in the sensible heat regenerator. It was first introduced by Banks [3] and Maclaine cross [4] then Jurinak [5] and Stabat [6] improved the model. The following assumptions are considered:

• The state properties of the air streams are spatially uniform at the desiccant wheel inlet

• The interstices of the porous medium are straight and parallel

• No leakage or carry-over of streams

• The interstitial air velocity and pressure are constant

• Heat and mass transfer between air and porous desiccant matrix is considered using lumped transfer coefficients

• Diffusion and dispersion in the fluid flow direction are neglected

• No radial variation of the fluid or matrix states

• The sorption isotherm does not represent a hysteresis

• Air reaches equilibrium with the porous medium

image667

The heat and mass conservation equations:

image668 image669 image670

Heat and mass transfer equations:

image671

Equations (2), (3), (4) and (5) are coupled hyperbolic non-linear. With the assumption of the Lewis number (Le), equal unity and the desiccant matrix in equilibrium with air means (Td= Teq and weq= wa). Banks [2] used matrix algebra and proved that these equations can be reduced using potential function Fi(T, w) to the following system:

(7)

9(C*)

(10)

Where

C * = Mm. c pm N

Подпись:Подпись: = 0

image674

Подпись: cПодпись: pm

image677

Introducing the equations of heat transfer alone in the sensible regenerator as stated in [3]:

(12)

ncf is the efficiency of the counter flow heat exchanger for balanced flow.

F = h

(13)

Подпись: F = (273.15 + T) F 2 =

image679

6360 (14)

The model described above, was implemented to SPARK [8] a general simulation environment based on equation. In the next section the model is validated experimentally.

Solar fraction and electricity consumption

From the 282 operation hours in the cooling mode about 35.8% (101 hours) no solar operation was possible. On 148 hours (52%) a complete solar operation covering the whole hour was possible, on the remaining 33 hours the solar system could only provide a part of the driving heat. From an energetic point of view, 59.8% of the driving energy over the whole cooling operation period came from the solar system.

The electricity consumption over the whole cooling period was 10% of the cooling energy. This electricity consumption includes the electricity needs for the chiller itself, the heat rejection via the boreholes, the driving circuit and the cold distribution to the cooling coil in the air handling unit. But it does not include the pump of the solar loop.

image632

temperature lift: T_MT_in — T_NT_out [K]

Подпись: driving temperature (T_HT_in) [°C]

Fig. 5. Frequency diagram of driving temperatures Fig. 6. Frequency diagram of temperature lifts.

(from heating net or solar system).

Calculating an electric COP for the whole system as the cooling energy produced per kWh of electric energy consumed, an electric COP of approximately 10 can be derived. It has to be noted however, that although energy efficient pumps have been used, all pumps operate at constant volume flows independently of the chillers output power. Thus a constant electricity consumption of about 500W is measured during operation. Optimising the control, e. g. using a flow rate control when low cooling capacity is needed may give potential for reducing the electricity consumption.

Experimental validation

image680

The experimental installation presented in section 2 is used to validate the above presented model. Different inlet conditions were considered with temperature varying form 25°C to 38°C and humidity ratio from 10 to 15 g/kg and for different regeneration temperature (55°C, 60°C 65°C, 70°C and 80°C).

The figure 2 [9] plots the temperature error versus the humidity ratio error for the most significant points. It shows the maximum predicted error 0.4 g/kg and root mean squared error (RMSE) is 0.22 g/kg. For the temperature and the maximum committed error in the outlet temperature is 1 °C for RMSE of 0.65 °C. It must be noted that the uncertainty in the measurement of the mean values is 0.7°C for the temperature and 0.3 g/kg. It shows that the committed errors are in the domain of the uncertainty of the measurements.

In the following section the experiments are performed on the desiccant using experimental measurements and numerical results in order to perform the sensitivity analsysis.

Description of the novel generator design

In this heat and mass transfer prototype, the desiccant solution and the air stream is brought into contact with each other in a cross flow configuration. The regenerator consists of a stack of polycarbonate (PC) twin wall plates. Each plate has an internal heating water circuit. The plates are covered with a wick fibre to facilitate intimate contact between the liquid desiccant and air. This is done to increase the exposure time over the plates and thereby enhance the desired mass transfer and heat exchange. The general design is shown in the left part of Fig. 2.

The advancement of the distribution system consists of:

• a separation between the liquid desiccant distributor and the contact area between air and desiccant. The LiCl solution is distributed over the textile separately in an advance stage i. e. before coming in contact with the air stream. A separator is used to split the generator into two chambers. In the upper small chamber LiCl throttled through perforated plexiglass tubes, completely diffused through the textile fibres and trickles down to come in contact with the horizontal air stream in the lower large chamber. This separator will prevent LiCl particles to drift into the air stream, resulting in a zero carryover regenerator.

• the employment of promising fibres, a new natural man made fibre produced from wood pulp. This fiber can absorb 50% more moisture than cotton did, it has a lower sorption index than cotton, which means higher transport speed of the liquid.

• the employment of perforated plastic tubes with different tested throttling-points diameters.

The LD distributor in this design uses a plurality of parallel plexiglass tubes to horizontally distribute the LD over the wick. The parallel tubes extend outwardly from openings in the lower edge of one of the sides of a LD feed box and are closed from the free end. The tubes penetrate the PC plates horizontally and deliver the desiccant solution over the coated plates at a number of equally spaced-apart locations (discharge-holes). The distributor is shown on the right part of Fig.

2.

The discharge holes are preferably formed on the plexiglass tubes by using a CNC machine to make the fine drip points with high precision. The size and number of the discharge holes are selected to provide the desired liquid flow. The throttling-point density determines the flow into the wick material covering the PC twin-wall plates. Likewise, the distance between the discharge — holes is selected to accommodate the desired LD flow rate with the maximum even distribution.

image618 image619
image620

The regenerator consists of an upper and a lower hot water-feeder box. Each water feeder box consists of four equal chambers separated by baffles. Hot water will pass through the internal passages of the PC twin-wall plate in a 3-loop serpentine path. It will enter the first chamber from a primary opening in the upper wall of the chamber, and it collides with a liquid spreading surface that faces the internal PC passages. The spreading surface will also disrupt the downward flow of water and will direct it equally to the channels.

Fig. 2. Exploded view of generator design along with the illustration of flow directions of heating water, LD and air (upper left), an isometric enlargement of the LD distribution system and a top-view of how the glass tubes penetrate the polycarbonate twin-plates (upper right).

Operation day in the cooling mode

In Erro! A origem da referenda nao foi encontrada. an example of the operation during a normal summer day is shown. The diagram corresponds to June 13th, 2007. Cooling operation starts at 7:30 and ends at 14:30. In the morning till 11:30 the chiller is operated with heat from the heating net as the temperature in the solar storage was not sufficient to operate the machine. The shaded area between 11:30 and 14:30 shows the time when the system was operated with heat from the solar system.

Large temperature variations can be observed in the driving circuit. While the cyclic temperature peaks observed in the return flow from the chiller in all circuits (T_HT_out, T_MT_out and T_NT_out) are characteristic of the periodic working adsorption machine, the variations in the driving inlet temperature (T_HT_in) are due to the dynamic behaviour of the plate heat exchanger within this circuit but also variations in the temperature of the heating net. These variations were not expected. A much more stable inlet temperature is observed during the solar operation period. Further it can be seen, that the boreholes effectively attenuate the temperature peaks giving smooth and constant feed temperatures (T_MT_in) to the machine.

4.1. Heating operation

For the heating operation the adsorption chiller is operated as a heat pump. For this operation the evaporator is connected to the borehole system and the heat rejection to the heating coil in the air handling unit.

Sensitivity results and analysis

4.1 Results

Experiments have been conducted on the desiccant wheel with the parameters (e. g. outside temperature, outside humidity ratio, regeneration temperature and regeneration humidity ratio) varying in the range defined in the section 2. The complete DOE of 4 parameters operating between 2 levels needs 24=16 experiments. When a combination of the studied parameters was

difficult to achieve experimentally (3 experiments) the results of the model was used to complete the DOE. Table 1 below shows the results

Table 1. Dehumidification rate of the desiccant wheel for different operating conditions

Ti

w1

T8

W8

w1-w2

T1

w1

T8

w8

w1-w2

[°C]

[g/kg]

[°C]

[g/kg]

[g/kg]

[°C]

[g/kg]

[°C]

[g/kg]

[g/kg]

25

11

55

10

4.8

35

11

55

10

2.9

25

11

55

15

3.6

35

11

55

15

1.7*

25

11

75

10

6.8

35

11

75

10

5.2

25

11

75

15

5.8

35

11

75

15

4.1

25

14.5

55

10

5.4*

35

14.5

55

10

4.1

25

14.5

55

15

4.3

35

14.5

55

15

2.95

25

14.5

75

10

8

35

14.5

75

10

6.7

25

14.5

75

15

7

35

14.5

75

15

5.6*

*

Calculated by

the model

image681
Solving 16 equations for the effects identified in the equation (1). The dehumidification of the desiccant wheel is then written in function of the operating parameters

4.2 Analysis

Form the equation identified we notice a constant dehumidifying effect of 4,987 g/kg which will be increased or decreased depending on the operating conditions. In order to compare the effect of each parameter with the constant effect, they are divided by this later and a graphical comparison can thus be established. The figure below shows the effect of each parameter.

As commonly known, figure 3 shows tha the regeneration temperature has the most significant impact on the dehumidification performance of the desiccant wheel. In the same time outside conditions has an important impact on too.

• Increasing the regeneration temperature from the mean value (65°C) to its upper limit will increase the dehumidification 23%

• Increasing the outside temperature to its higher limit will decrease the dehumidification of 17%.

• Increasing the outside humidity ratio will increase the performance of the desiccant wheel of 13%

• Increasing the regeneration humidity ratio will decrease the performance of the desiccant wheel of 8%.

The reason behind these observations is the vapor pressure difference of the air and that of the surface of the silica gel. This vapor pressure difference is the driving force of the adsorption phenomena. Reminding that if the desiccant is cold the vapor pressure at its surface. For the moist air the temperature does not have a significant impact on the vapor pressure while humidity ratio does.

The desiccant coming for the regeneration sector is hot and dry, it is first cooled down by the outside air in the process sector and then the adsorption occurs.

image682

Fig. 3 Effect of the operating conditions and their combinations

So if we increase the temperature of the outside air its vapor pressure will not increase significantly in reversal it will not cool sufficiently the hot desiccant form arriving from the regeneration sector yielding a high vapor pressure at the surface of the desiccant and thus the adsorption will be less efficient. If we increase the humidity ratio of the outside air its vapor pressure will increase significantly yielding a higher vapor pressure difference thus increasing the adsorption. When increasing the regeneration temperature, the regeneration air vapor pressure does not increase significantly while the desiccant is heated and its vapor pressure increases leading to an efficient drying of the desiccant. This desiccant leaving the regeneration sector is dry and thus having fewer vapor particles, yielding a low vapor pressure at its surface and thus high adsorption. When increasing the regeneration humidity ratio we will increase the vapor pressure of regeneration air stream, reducing the transfer from the desiccant to the regeneration air, yielding less drying of the desiccant. When the sorbent will leave the regeneration with more water vapor at its surface its vapor pressure is then higher, reducing thus the adsorption.

This clearly shows the limitations of the desiccant cooling technique regarding outside conditions. It demonstrates that high outside temperature reduces significantly the performance of the desiccant wheel. Regarding the outside humidity ratio even if the dehumidification increase with increasing outside humidity ratio, we noticed that for outside temperature beyond 30°C the maximum dehumidification rate is 6 g/kg. Taking into account the maximum humidity inside the building (e. g. 11.8 g/kg) and the humidification across the supply humidifier we conclude that the maximum outside humidity under which a desiccant system will operate efficiently is 14.5 g/kg.

3. Conclusion

In this paper a sensitivity analysis based on the design of experiments was conducted on a desiccant wheel using mainly experimental results. A numerical model was validated experimentally and provided the missing combination of the experiments. The impact of the outside and regeneration conditions on the dehumidification rate of the desiccant wheel was studied. As widely known the regeneration temperature has the most significant impact on the dehumidification rate, but the impact of the outside temperature, outside humidity ratio and the

regeneration humidity ratio are very important as well. These results showed that desiccant cooling is an interesting option for moderately hot and moderately humid climates

Nomenclature

Cpa

heat capacity of air [J. Kg-1.K-1]

NUT

number of transfer unit [-]

cpv

heat capacity of water vapour

RMSE

Root mean squared error

[J. Kg’.K-1]

[unity of the variable]

cpm

heat capacity of the regenerator matrix

t

time [s]

[J. Kg-‘.K-1]

T

temperature [K]

C

heat capacity (J. Kg"1.K"1)

T

a

air temperature [K]

Fi

potential characteristic [-]

T

eq

equilibrium temperature [K]

ha

enthalpy of moist air [J. Kg-1]

Tm

matrix temperature [K]

H

enthalpy of the desiccant [J. Kg-1]

u

fluid velocity (m. s-1)

Jt

lumped heat transfer coefficient [s-1]

Wa

humidity ratio of moist air. [Kg/Kg]

Jm

lumped mass transfer coefficient [s-1]

Wd

water content of desiccant [Kg/Kg]

Le

Lewis number [-]

z

coordinate in the flow direction [m]

ma

air mass flow rate [Kg. s-1]

T

ratio of matrix mass over air mass [-

Md

mass of the desiccant matrix [Kg]

n

efficiency [-]

Mm

mass of the aluminium matrix [Kg]

p

density [Kg. m-3]

N

angular speed of the wheel [rd. s-1]

Yi

parameter [-]

Acknowledgments

The authors would like to thank Mr Michel Burlot for his valuable technical support on the experimental

installation. This work was supported by ADEME (French Agency of the Environment and Energy

Management) and the regional council of Poitou Charentes.

References

[1] P. Bourdoukan., E. Wurtz., P. Joubert., M. Sperandio,. (2008). Potential of solar heat pipe vacuum collectors in the desiccant cooling process: modelling and experimental results. Solar Energy. 0.1016/j. solener.2008.06.003

[2] J. Goupy, (2001). Introduction aux plans d’experiences, Dunod.

[3] P. J. Banks, (1972). Coupled equilibrium heat and single adsorbate transfer in field flow through a porous medium — 1. Caracteristics potentials and specific capacity ratios. Chemical Engineering Science, 27, 1143-1155.

[4] I. L. Maclaine-cross, P. J. Banks, (1972). Coupled heat and mass transfer in regenerators — prediction using analogy with heat transfer. Int. J. Heat Mass Transfer, 15, 1225-1242.

[5] J. J. Jurinak. (1982). Open cycle solid dessicant cooling — component models and system simulation, Ph. D, Winconsin, Madison.

[6] P. Stabat, (2003). Modelisation de composants de systemes de climatisation mettant en oeuvre l’adsorption et l’evaporation d’eau, Ph. D, Ecole des Mines, Paris.

[7] W. M. Kays, A. L London, (1984). Compact Heat Exchangers, McGraw-Hill, New York.

[8] SPARK. (2003). Simulation Problem Analysis and Research Kernel. LBNL, Berkeley, California.

[9] P. Bourdoukan., E. Wurtz., P. Joubert., M. Sperandio, (2008). Critical efficiencies of components in desiccant cooling cycles and a comparison between the conventional mode and the recirculation mode."

In proceedings of ECOS conference, Krakow, Poland, 1, 435-443

Experimental setup

Two different sets of experiments were carried out to prove the advancements on the distribution system.

Подпись: Fig. 3. Experimental setup to determine the optimum size and number of discharge points for the distributor (upper), and to check the LiCl diffusion through the textile (lower).

1. The objective of the first test is to determine the optimal size and number of discharge points distributed along the plastic tube. The optimized tube will serve the LD flow rate required by the application and ensure fairly equal amounts of the LD throttled-out through each discharge bore. Five perforated plastic tubes are tested. Each tube is perforated by a CNC machine in equally-spaced positions. H2O/LiCl (30% salt concentration by weight) is used as the liquid desiccant with three volume flow rate values. In order to investigate the effect of the distance between the LiCl feed-box and the discharge-point location, different positions along the tube are selected. Transparent PVC hoses are used to enclose the discharge-bores, without affecting the throttled LiCl. The free ends of the PVC hoses are connected to plastic bottles in order to collect the LiCl throttled out of the intended discharge-bore. Fig. 3 shows the experimental setup.

2. The second test is a comparison between different types of textiles (100% cotton, 100% viscose, 100% polyester, 100% polyamide and wood pulp based textiles) with different thicknesses and different weaving. Each type of the mentioned textiles has been tested to measure the absorption capacity and the diffusion speed by simply pouring a quantity of LiCl onto the upper surface of a taut piece of textile, and measuring the needed time for the LiCl droplets to be completely absorbed by the textile fibers. Each type has been tested in both, dry and wet state. Then PC plates were coated with the textiles which show the best absorption and diffusion speed, and exposed to a LiCl solution throttled through the perforated plexiglass tube. A violet fluorescent light and dyes have been used to support the visual inspection of LiCl distribution through the textile fibers. Fig. 3.

Heating power and COP

As for the cooling operation, for the following statistics on hourly values, only hours with a constant operation of the machine were considered in order to be sure to consider steady state conditions only. This covers more than 63% of the hours with operation, as can be deduced from Table 1. Fig. 8 shows the frequency diagrams of the hourly COP and heating powers (Fig. 9). For the whole heating operation period a thermal COP of 1.43 is obtained. The hourly mean heating power was 9.43kW with a standard deviation of 1.59kW.

image634

in the period between 11:30 and 14:30.

Simulation Model and Validation

A detailed full dynamic simulation model of the installed system which considers the electricity consumption of all installed components (fans, pumps, etc.) has been developed by zafh. net in the simulation environment INSEL [5]. The simulation model is used to analyse the effect of different storage charge and discharge options. A very small time step of 10 seconds is used for an accurate consideration of all thermal capacities in the complex system. The complete model has been validated against measured data of the installed system. Figure 2 shows the measured performance of the solar cooling system together with the simulation results for one day in August 2007 [1]. A comparison of the simulated and measured outlet temperatures of the generator, condenser and evaporator of the ACM and of the collector field clearly visualises that the performance of the installed system is very well described by the developed simulation model.

[1] Thermo-economic optimisation shows that for solar thermal systems it is beneficial to oversize the solar

field.

[3]

The design and method of manufacture has been protected internationally by patent applications which have been lodged in Patent Cooperation Treaty member countries.

[4] Introduction

Solid desiccant cooling has been proposed as an alternative to vapour compression refrigeration for space cooling. It is an environmentally attractive solution, which does not require ozone depleting refrigerants and can be run off low temperature waste heat or solar heat.

Utilising a solar heat source appears to offer potential for approaching zero emissions HVAC. However, some form of thermal back up, such as an electric chiller or a gas heater, is normally provided to overcome the intermittent availability of solar heat. Significant fossil fuel energy consumption and associated greenhouse gas emissions can arise from this backup thermal source.

To minimise greenhouse gas emissions, it is useful to explore the potential for operating without backup, and allowing increased fluctuations in thermal comfort conditions inside the occupied space of the building. In a previous paper, White et. al. [1] investigated the performance of a once — through solar desiccant cooling system, for airconditioning a commercial office space. The study utilised the TRNSYS computer simulation software package to explore the range of thermal comfort conditions that might be expected inside a typical office building when the solar desiccant system is operated, without a backup supply of cooling or heating, in Sydney, Australia.

Using an airflow optimised for a 60°C regeneration air temperature (4.53 airchanges per hour), the temperature inside the occupied space was found to exceed the target temperature of 26°C for 196, 131 and 87 hours per year for flat plate solar collector areas of 0.083, 0.167 and 0.333 m2 per m2 of floor area respectively. These results appear to indicate that some form of backup thermal source is likely to be required for a small number of hours per year.

Analysis of the results

The first set of experimental data for the distribution tube was analyzed by the factorial design analysis using MINITAB version 15.

Factorial design analysis is a statistical method in which every level of one factor is tested in combination with every level of another factor. In general, in a factorial analysis, all possible combinations of factor levels are tested [5].

Four variables (factors) were considered to have the highest influence on the distribution test. The first factor is the LiCl volume flow rate with three treatments (levels) with a range between 0.3 l/min-1.2 l/min. The second factor is the bore size with five levels, ranges between 0.5 mm — 0.9 mm. The third factor is the distance for each throttling-bore from the LD feed-box entrance with ten levels, ranges between 4 cm-58 cm. The forth factor is the run time for each test with two levels (10 min and 15 min). Two replicates were taken for each trial in order to get an unbiased estimation. Since the number of levels for the studied factors is not equal, general full factorial design was used to analyze the measured data. The significance level a of the test is 0.05.

Подпись: Fig.4 depicts the normal probability plot of LiCl mass throttled through the perforated plastic tubes. The normal probability plot is a graphical technique for assessing whether or not a data set is approximately normally distributed. The data are plotted against a theoretical normal distribution in such a way that the points should form an approximate straight line. Departures from this straight line indicate departures from normality. The figure shows a little bit lower tail that do not fall exactly along a straight line passing through the centre of the plot. This indicates some potential problem with the normality assumption. However, no severe deviations from normality are obvious apparent. It is clearly shown that the minimum mass flow rate that could be used to ensure normal probability assumption corresponds to a bottle content with a total mass of 230 g.

From the output a large F-statistic and a low p-value is seen, indicating strong evidence of a correlation between the four factors and the LiCl throttled through each discharge-point. It is also shown that p-values for all terms in the test are less than the significance level. This confirms the quality of the test.

Fig. 4. Normal probability plot of the response (LiCl mass)

Fig. 5. Residual versus the response (LiCl mass)

Подпись: Fig. 6. LiCl adsorption-speed test for different textiles

The second set of tests focused on the absorption capacity and diffusivity of different textiles. A probe of lithium chloride solution (45% salt concentration by weight) was poured on a taut piece of textile. Then the needed time for the LiCl droplets to be completely adsorbed by the textile fibers were measured. The tests were carried out for the textiles in dry and wet state. Fig. 6 shows the adsorption speed for probes of one of the promising wood-pulp based textiles, polyamide, viscose, cotton, polyester and polystyrene fibers.

The best results were obtained for pure wood-pulp, polyamide and polyester fibres. Different probes of polyester and polyamide fibres reveal essentially inferior adsorption speed. Therefore the material itself is not the only factor influencing the adsorption speed. During the tests a strong dependency was detected between the sorption speed and the material thickness, the pore size and the weaving structure. Further tests must be carried out to clarify these influences. A good LiCl diffusion behaviour through the pure wood-pulp fiber was proved in the experiments.

2. Conclusion and Outlook

This paper presents experimental tests and a factorial design analysis of a desiccant distributor system which will be used in the prototype of a LD regenerator. It shows the optimal conditions that will ensure even distribution and the applicability to serve a wide range of desiccant flow rates. Furthermore, it presents experimental investigations for different textiles in order to choose the best textile with the highest adsorption speed and the finest desiccant diffusion through the textile fibers.

The factorial design analysis shows the optimal diameter of the throttling-holes, and the optimal spacing between those throttling diameters. Those optimal values will ensure a fairly even distribution and it will serve a wide range of desiccant volume flow rates, which will make the LD distributor flexible to serve the requested volume flow rate requested by the air conditioning cooling load.

The experimental investigations related to textiles approved that new textiles such as pure wood — pulp fibers and polyamide perform much better than traditional textiles currently in use such as cotton.

Further investigations must be carried out to clarify the influence of the textile thickness, the pore size, the desiccant type and the weaving structure on the adsorption behaviour of textiles. Experimental investigations might be considered concerning: the durability of the throttling-bores against the desiccants corrosivity, the influence of the regenerator operating conditions (high temperatures and corrosivity) on the textile performance, the effect of using higher desiccant concentration e. g. 40% on the distribution through the perforated tubes (in case if the LD distributor will be used in a conditioner instead of a regenerator), the necessity to flush the distributor and the solution that might be used.

References

[1] Kathabar technical information, www. kathabar. com.

[2] Krause M., Saman W., Vajen K. (2005), Regenerator Design for open Cycle Liquid Desiccant Systems — The-oretical and Experimental In-vestigations, Proc. Interna-tional Conference Solar Air-Conditioning, Staffelstein, 6. — 7.10.2005.

[3] Lowenstein A., Slayzak S., Kozubal E. (2006), A Zero Carryover Liquid Desiccant Air Conditioner for Solar Applications, ASME/Solar06, Denver, USA.

[4] Lavemann E., Peltzer M. (2005), Solar Air Condition-ing of an Office Building in Singapore using Open Cycle Liquid Desiccant Techno-logy, Proceeding of the International Conference on Solar Air Condi-tioning, Staffelstein, 06.-07.10.2005

[5] Montgomery, D., Runger G. (1999), Applied statistics and probability for engineers, second edition, John Wiley and Sons, Inc. New York.