Modeling of the Pyrolysis of Large Samples of Polyethylene Including

The polyethylene undergoes important volume variations during the experiments: expansion during the melting process and contraction due to pyrolysis...
0 downloads 0 Views 141KB Size
436

Energy & Fuels 2002, 16, 436-442

Modeling of the Pyrolysis of Large Samples of Polyethylene Including the Melting Process J. Ceamanos,* J. F. Mastral, and F. Liesa Department of Chemical and Environmental Engineering, University of Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza, Spain Received July 17, 2001. Revised Manuscript Received November 20, 2001

An experimental and theoretical study has been carried out on the thermal degradation of large samples of high density polyethylene in an inert atmosphere. The experimental results show the existence of significant longitudinal temperature profiles inside the solid, which are basically due to its low thermal conductivity and the melting process. The polyethylene undergoes important volume variations during the experiments: expansion during the melting process and contraction due to pyrolysis. The mathematical model developed takes into account the volume variation and the melting process using the enthalpy method. The temperature profiles and the solid conversion are predicted, and good agreement with the experimental temperature results is obtained. Some discrepancies are observed in the case of solid conversion, probably due to the resistance to mass transfer of product gases in the solid.

Introduction The theoretical study of the physicochemical processes occurring during the heating of a solid undergoing an incident heat flux is of considerable interest for different applications such as the knowledge of the minimum residence time necessary for the calculation of the global conversion of the solid or for the volatile emission. The solid conversion modeling allows the estimation of the time necessary for a given conversion and the prediction of the volatile emission rate. This prediction, combined with the study of the degradation of the material in an air atmosphere can be useful for the calculation of the time to ignition or the combustion rate of a material in the event of a fire.1,2 The volatile emission rate obtained from the modeling of the pyrolysis of the material should also be included in the study, since in some cases the combustion rate is controlled by the volatile supply from the solid thermal decomposition. The surface of the material undergoes combustion but when the access of oxygen to the interior of the sample is restricted, it degrades under pyrolytic conditions. In a previous work3 a model was presented for the prediction of temperature profiles and solid conversion of a material exposed to a heat flux in inert and oxidative atmospheres. The model was experimentally validated for different types of wood and checked with other materials such as polyurethane foam. In the * Corresponding author. Fax: (+34) 976 762142. E-mail: ceamanos@ posta.unizar.es. (1) Bilbao, R.; Mastral, J. F.; Lana, J. A.; Ceamanos, J.; Aldea, M. E.; Betra´n, M. A model for the prediction of the thermal degradation and ignition of wood under constant and variable heat flux. J. Anal. Appl. Pyrol. 2001, 62 (1), 63. (2) Bilbao, R.; Mastral, J. F.; Lana, J. A.; Aldea, M. E.; Ceamanos, J.; Betra´n, M. Experimental and theoretical study of the ignition and smouldering of wood including convective effects. Combust. Flame 2001, 126 (1-2), 1363. (3) Bilbao, R.; Mastral, J. F.; Ceamanos, J.; Aldea, M. E. Modeling of the pyrolysis of wet wood. J. Anal. Appl. Pyrol. 1996, 36, 81.

comparison of the predictions with the results of the experiments performed, the structural stability of the material was shown to be very important in the obtention of a suitable agreement. For wood, which basically maintains its structure during its thermal decomposition, a good agreement is observed. However, when applied to materials experiencing expansion (due to density differences during melting) and contraction of the physical structure (due to thermal decomposition), the model, which does not include dimensional variations, cannot be considered valid. When studying materials such as polyolefines and other polymers, other important aspects influencing the model results that need to be addressed are the melting over a temperature range and not a specific temperature, and the variation of the thermal properties from one phase to the other. In general, a semicrystalline thermoplastic such as HDPE can be described as a mixture of amorphous and crystalline fractions. The melting point is defined as the temperature at which an isolated crystal loses its crystallinity.4,5 This temperature equals 414.6 K for isolated crystals. However, this value can be affected by different causes. The melting temperature has been measured6 for different types of polyethylene crystals undergoing different heating rates. A decrease of the melting point of around 8 °C was observed for heating rates of 15 °C/min, although no significant changes were observed for higher heating rates. In the case of a large sample of the material, there are no isolated crystals, but a mixture of the above-mentioned fractions, and numerous amorphous-crystalline subsystems with dif(4) Wunderlich, B. The Melting of Defect Polymer Crystals. Polymer 1964, 5, 125, 611. (5) Wunderlich, B. The basis of thermal analysis. In Thermal Characterization of Polymeric Materials; Turu, A., Ed.; Academic Press Ltd.: San Diego, 1997; Vol. 1. (6) Hellmuth, E.; Wunderlich, B. Superheating of linear highpolymer polyethylene crystals. J. Appl. Phys. 1965, 36 (10), 3039.

10.1021/ef010179e CCC: $22.00 © 2002 American Chemical Society Published on Web 02/05/2002

Melting and Pyrolysis of Polyethylene

Energy & Fuels, Vol. 16, No. 2, 2002 437

Figure 1. Experimental system.

ferent energetic states. In this case, melting consists of a decrease of the crystalline fraction from its original value in the sample down to zero. Moreover, when a relatively large sample is heated, different heating rates are usually observed at different regions of the sample. Taking this into account, the factors affecting the initial crystalline fraction, the characteristics of the crystals, and the heating rate would all affect melting, and this will occur over a temperature range.5 Thus, a study of melting is necessary in which a range of temperatures is considered. Taking this into account, the objective of the present work is the modeling of the melting of relatively large samples of high-density polyethylene, its influence on the heat transfer, and the thermal decomposition of the material. A mathematical model of the process is presented in which the material undergoes different heating rates varying over time. It predicts the duration of melting, the weight loss and the variation of the sample thickness (expansion and to shrinking). Experimental Section The experimental system is shown in Figure 1. The sample is heated by a 200 × 200 mm square radiant panel inside a 600 × 500 × 560 mm refractory steel chamber. The sample holder and the radiant panel are arranged in parallel and horizontally. Several thermocouples are placed in the sample at different depths from the base surface. The experimental system also allows continuous measurement of the solid weight by means of a precision balance connected to the sample holder. The samples consist of 50 × 50 × 20 mm slabs of highdensity polyethylene. The thermal decomposition of the sample

is carried out in an inert atmosphere. A flow rate of 1.2 m3 h-1 of nitrogen was used. Only one side of the sample is exposed to the radiant heat flux. This arrangement provides a uniform heat flux over the sample surface so that heat transfer is one-dimensional inside the sample. To corroborate this behavior, some experiments were carried out with thermocouples placed at different points of the sample surface. No significant differences between the temperatures measured at different points was found. The experimental system may show some similarities with the Cone Calorimeter (ASTM E-1354-99, ISO5660), widely used to measure heat release rate, time to ignition. The dimensions of the sample are similar and the heat flux over the sample is basically radiant in both experimental systems. However, there are some differences which should be noted. Cone calorimeter experiments are carried out in an air atmosphere. In the system used in this work, experiments under different atmospheres (inert, air, and atmospheres with different oxygen concentrations) can be carried out. Cone calorimeter experiments are carried out under constant radiant heat fluxes, while a variable heat flux has been used in the system used in this work. Variable heat flux can be considered more realistic in the event of a fire, but presents more difficulties when trying to describe the process theoretically. Experiments for four different radiant heat fluxes, qe′′ were carried out. These fluxes were experimentally measured using a Gardon-type, water-cooled and gas-purged heat flux gauge and the curves obtained are shown in Figure 2. As can be observed, the heat fluxes are variable over time and, for a given power supply, increase up to a constant value. These maximum values range between 5 and 35 kW/m2. Figures 3 and 4 show the variation during the experiment of temperatures measured at 1, 5, 10, 15, and 20 mm from the sample bottom for heat flux curves 1 and 2. As will be

438

Energy & Fuels, Vol. 16, No. 2, 2002

Figure 2. Radiant heat fluxes vs time.

Figure 3. Experimental temperature profiles vs time for heat flux curve 1.

Figure 4. Experimental temperature profiles vs time for heat flux curve 2. remarked later on, given that the surface position changes over time, the location of the different measuring points in the sample thickness are considered from the base of the sample. As can be observed, the temperature increases up to a value between 400 and 415 K when a plateau appears that extends for longer as the distance from the sample surface increases. The heat transferred to the region of the sample where melting occurs is consumed in this process, so no temperature increase is observed. The melting process increases the temperature gradients throughout the sample thickness. The thermal conductivity of melted polyethylene is lower than that of solid,

Ceamanos et al.

Figure 5. Experimental and predicted global solid conversions vs time for heat flux curves 2, 3, and 4. so the temperature gradients in the sample are higher. Another important cause of the significant increase in the temperature gradients is that the variable incident heat flux over the surface increases during the experiment. The melting effects on the heat transfer behavior are more evident for low heat fluxes, as can be clearly observed in Figure 3. For a slghtly higher heat flux, as the one shown in Figure 4, it can be observed that pyrolysis of the surface starts approximately when melting of the bottom surface has finished. For higher heat fluxes the effects associated with melting can be hidden by the weight loss and volume reduction of the sample due to thermal decomposition. For these heat fluxes the pyrolysis can start at the surface when part of the sample is still undergoing melting. Once the pyrolysis temperature is reached, the material generates basically gases and vapors. A regression of the surface starts so that the heat flux incides more directly onto inner regions of the sample, increasing the rate of the global process, accelerated by a decreasing mass of remaining material. The pyrolysis is endothermic but its enthalpy is relatively low compared to the radiant heat flux, so this heat consumption affects the temperature profiles only slightly. Figure 5 shows the variation of solid conversion during the experiment for heat fluxes 2, 3, and 4. The theoretical conversions predicted with the model are also shown, although they will be discussed below. The conversion for heat flux 1 is not shown since the maximum value observed is too low (lower than 10%), and is obtained over a very long time (2.5 h). Obviously the time necessary for the onset of the pyrolysis decreases as the heat flux increases. Once the surface reaches the temperature of pyrolysis, this takes place over a narrow range of temperatures. The rate is higher for heat fluxes 3 and 4 than for 2, but is similar for 3 and 4. This could indicate that the pyrolysis rate is controlled by a phenomenon other than heat transfer for high heating rates.

Mathematical Model The mathematical model used describes the thermal decomposition of a material exposed to a variable heat flux. It includes the basic physical and chemical phenomena taking place during the thermal decomposition of the material and does not use fitting parameters. The main assumptions of the model are the following: • The solid surface is exposed to a variable radiant heat flux. Heat losses due to re-radiation and conduction toward the solid are taken into account. Convective losses and the effect of the evolved gases on the heat transfer at the surface are considered negligible. • Heat transfer in the solid is one-dimensional. The dependencies of thermal conductivity and the heat capacity on temperature are taken into account.

Melting and Pyrolysis of Polyethylene

Energy & Fuels, Vol. 16, No. 2, 2002 439

• No resistance to mass transfer of the gases formed during the thermal decomposition is considered. Gases leave the solid immediately after they are formed. • The kinetics of pyrolysis, obtained by thermogravimetric techniques, is considered to follow a first-order power law. • The melting of the polyethylene is considered to occur over a range of temperatures. The dependency of melting enthalpy on the temperature is included. During melting the crystalline fraction depends lineally on temperature. • The volume of the sample changes due to both melting and pyrolysis. The surface regression is predicted by the model. This transient system can be described by the following time-dependent partial differential equation for onedimensional problems:

∂(FH) ∂T ∂ K + (-∆Hr)(-ra) ) ∂t ∂x ∂x

(

)

(1)

where F is the density of the material, H is the enthalpy of the material, K is the thermal conductivity, T is the temperature, (-∆Hr) is the enthalpy of thermal decomposition, and (-ra) is the reaction rate. The left-hand side, the variation of density and enthalpy, takes into account the variation of sensible heat, density, and the melting enthalpy. The second term on the right-hand side, the heat consumed by the thermal decomposition, becomes zero for temperatures lower than that of pyrolysis. The mass balance describes the weight loss due to the thermal decomposition:

(-ra) ) -

dXs ∂m ) m0 ∂t dt

(2)

where Xs is the solid conversion, defined as Xs ) (mo - m)/mo. The boundary conditions are as follows:

t ) 0, ∀x

T(x) ) T0 Xs ) 0

(3)

t > 0, x ) 0

∂T )0 ∂x

(4)

t > 0, x ) L

qi′′ - σT4S - K

∂T )0 ∂x

(5)

Equations 4 and 5 describe the heat balances for the lower and upper surfaces of the sample, respectively. The incident radiant heat flux was determined previously in the preliminary experiments, and the radiant heat losses were calculated using an emissivity of 0.5 for the material and the surface temperature calculated from the model. Model parameters: The parameters used to solve the equations have been obtained by using other experimental systems or from the literature. The melting temperature depends on the characteristics of the polymer crystals, the production process, and the heating rate of the material. In a relatively large sample, the heating rate can vary throughout the thickness, so the crystals closer to the sample surface undergo more severe heating conditions. However, this effect has not been considered and the same melting temperature interval is considered for the solid as a

whole. The loss of crystallinity and hence the melting is assumed to occur between 400 and 415 K, similar to the values observed for commercial HDPE.7 A sensitivity analysis has been carried out to determine the importance of this parameter and its influence on the predicted results. The polyethylene density decreases with increasing temperature, with a sharp increase during melting since the density of the crystalline phase is much higher than that of the amorphous phase. A density decrease obviously implies a volume increase, which was observed experimentally and estimated at approximately 16%. This significant variation of the property makes its description necessary in the model. The relationship between density and temperature takes into account the crystalline and amorphous phases, after the dependencies of the specific volumes on temperature:8

v ) xcvc + (1 - xc)va ) xc(0.993 + 3.0 × 10-4T) + (1 - xc)(1.152 + 8.8 × 10-4T) (6) where xc is the crystalline fraction of the material. For polyethylene at ambient temperature, a density value of 960 kg m-3 has been measured, so an initial crystalline fraction xc ) 0.66 is obtained. This relationship is valid only until melting, where the crystalline fraction decreases. The variation of the crystalline fraction, and hence the loss of crystallinity, is assumed to occur lineally over the range of temperatures where melting takes place, 400-415 K:

xc ) 0.66(415 - T)/15

(7)

The values of melting enthalpy and thermal conductivity also take into account the fraction of crystals present in the solid. Different authors9,10 have obtained the enthalpy of the material over a very wide range of temperatures, both experimentally and theoretically. The values of the thermal conductivity used in this work are Kcryst ) 0.45 W m-1 K-1 and Kamor ) 0.35 W m-1 K-1. The sample conductivity has been obtained from these values taking into account the crystalline fraction (x ) 0.66). During the range of melting temperatures, the value changes from that of the mixture of crystalline and amorphous fractions to that of the amorphous fraction, corresponding to a fraction of crystallinity equal to zero. The kinetic equation for the thermal decomposition of HDPE was determined previously11 assuming a simple first-order degradation:

dXs ) k(1 - Xs) dt

(8)

where k is the rate constant obtained in thermogravi(7) Ramos, M. A.; de Marı´a, M. R. Ingenierı´a de los materiales pla´ sticos; de Santos, D., Ed.; Dı´az de Santos: Madrid, 1992. (8) Chiang, R.; Flory, P. J. J. Am. Chem. Soc. 1961, 83, 2857. (9) ATHAS Databank, http://web.utk.edu/∼athas/databank/ intro.html 1990. (10) Gaur, U.; Wunderlich, B. Heat capacity and other thermodynamic properties of linear macromolecules. II. Polyethylene. J. Phys. Chem. Ref. Data 1981, 10 (1), 119. (11) Ceamanos, J.; Mastral, J. F.; Millera, Aldea, M. E. Kinetics of the pyrolysis of high-density polyethylene. Comparison of isothermal and dynamic experiments. J. Anal. Appl. Pyrol., in press.

440

Energy & Fuels, Vol. 16, No. 2, 2002

Ceamanos et al.

Table 1. Kinetic Parameters at Different Heating Rates T < 430 °C differential

T g 430 °C differential

β (°C min-1)

ko (min-1)

Ea (kJ mol-1)

ko (min-1)

Ea (kJ mol-1)

5 12 25 50

3.39 × 1018 9.34 × 1010 4.03 × 1010 5.08 × 1010

273.1 171.1 164.9 163.6

1.94 × 1026 1.14 × 1026 2.08 × 1025 1.55 × 1024

376.6 373.9 363.2 344.3

metric isothermal and dynamic experiments. The values of the rate constant previously determined are shown in Table 1. Model Solution: This problem can be described as a melting process with heat transfer in both the solid and the melted phases. When the material has a specific melting point, a defined surface between these two phases exists. Due to the calorific requirement for melting, the temperature gradient is not defined for this surface and eq 1 cannot be solved continuously for the whole sample. To obtain the solution of the heat diffusion equation in melting problems, two types of methods exist. The first type focuses on the obtention of the exact position of the interface by dividing the region into two subregions where the equation is solved by applying a moving boundary condition (Stefan condition) at the interface. The second type, fixed domain methods, solve the problem throught weak formulation,12,13 the so-called enthalpy method being the one used in this work. The phase change does not occur at a specific temperature, but over a range of temperatures, producing a so-called mushy region.14,15 This behavior makes the classic solution for the Stefan problem inappropiate, since the boundary conditions used imply a single phase change temperature. However, the enthalpy method allows the use of a range of melting temperatures, and it is relatively simple to implement. Different explicit and implicit numerical methods for the solution of the enthalpy formulation are available in the literature.14,16-19 Given the relatively short time required for the processes of heating, melting, and thermal decomposition, an explicit method has been used14 and has been shown to give adequate results. To obtain the solution of the heat diffusion equation, a variable change needs to be made:

Cp(T) dT ) dH

(9)

During melting the increase of the energy content with temperature must be added, taking into account (12) Oleinik, O. A. A method of solution of the general Stefan problem. Sov. Math. 1960, 1, 1350. (13) Friedman, A. The Stefan problem in several spaces variables. Trans. Am. Math. Soc. 1968, 133, 51. (14) Atthey, D. R. A Finite Difference Scheme for melting problems. J. Inst. Math. Appl. 1974, 13, 353. (15) O ¨ zis¸ ik, M. N. Heat Conduction; John Wiley & Sons: New York, 1980. (16) Meyer, G. H. The numerical solution of Stefan problems with front-tracking and smoothing methods. Appl. Math. Comput. 1978, 4, 283. (17) Voller, V.; Cross, M. Accurate solutions of moving boundary problems using the enthalpy method. J. Heat Mass Transfer 1981, 24, 545. (18) Whyte, R. E. An enthalpy formulation of the Stefan problem Siam J. Numer. Anal. 1982, 19, 1129. (19) Voller, V. Fast implicit finite-difference method for the analysis of phase change problems. Numer. Heat Transfer, Part B 1990, 17, 155.

Figure 6. Sensitivity analysis for the range of melting temperatures. Experimental and predicted temperatures considering melting temperatures 400-405 K, 400-410 K, 400415 K, and 400-430 K.

that the enthalpy of a finite volume is the sum of the enthalpies of both phases and the energy content due to melting. For the description of the volume variation, an invariable fixed grid with a given number of nodes is used in the calculations. Depending on the volume occupied by the sample, the boundary conditions will be applied to a different node, so that the sample surface position will change over time. This variation of volume is due both to the expansion process (which is taken into account with density variation) and to surface regression (due to pyrolysis which is taken into account through mass decrease). Sensitivity analysis: The melting process takes place over a range of temperatures which depends basically on the heating rate of the material. For higher heating rates, the ranges are wider. The same range has been used for the different heating rates used in this work. A sensitivity analysis has been carried out in order to determine the importance of this parameter on the predicted results. The ranges studied were 400-405 K, 400-410 K, 400-415 K, and 400-430 K. The temperature results show that only the temperature of the bottom of the sample is influenced by this parameter. In general, the increase of the temperature range gives better results as the heating rate increases. As an example, Figure 6 shows the variation of the temperature at this point for heating curve 3. Some differences are observed for the different conditions but only during the melting process. No differences are observed afterward. As can be observed, the optimum range for this heating rate would be between 15 and 30 K. For the heat flux curves 1, 2, and 4, the best ranges would be 15, 15, and 30 K, respectively. However, given the slight differences observed between experimental and theoretical results, a single range of melting temperatures (400-415 K) can be used. Discussion The model has been validated experimentally by comparing the results of temperature and solid conversion for the different heat flux curves used. Figures 7 and 8 show the comparison of the experimental and

Melting and Pyrolysis of Polyethylene

Figure 7. Experimental and theoretical temperature profiles vs time for heat flux curve 3.

Figure 8. Experimental and theoretical temperature profiles vs time for heat flux curve 4.

theoretical temperature profiles obtained for heat flux curves 3 and 4, respectively. For heating curve 3, Figure 7, a good agreement is observed. The duration of the plateau due to melting, observed especially for the bottom surface, is adequately predicted. The temperatures increase rapidly afterward up to a value close to 700 K. Figure 8 shows the comparison for the highest heating rate. As can be observed, a good agreement is also obtained. The time required for melting to occur is slightly overestimated by the model, as is indicated by the temperature of the bottom of the sample. This difference could be due to the adiabatic condition assumption considered in the model, which may not fit the experimental conditions. The temperatures increase rapidly after melting, obtaining very good agreement in general. Some differences are observed again at the bottom, where the predicted temperature increases very sharply. This increase starts approximately at 21 min and could be due to the small amount of remaining material. The experimental temperature increases more slowly probably due to the greater distance between the heater and the sample remains. In fact, this can be corroborated by plotting the thickness variation during the experiment, as shown in Figure 9, for the different heating rates used. The thickness increases from the initial value up to a maximum, subsequently decreasing to zero due to the thermal decomposition yielding

Energy & Fuels, Vol. 16, No. 2, 2002 441

Figure 9. Theoretical sample thickness vs time for heat flux curves 2, 3, and 4.

Figure 10. Theoretical sample thickness, surface, and bottom temperatures vs time for heat flux curve 2.

basically gases. When comparing the thickness variation for different heating rates, a decrease of the value of the maximum is observed when the heating rate increases. The maximum for the highest heating rate is sharper, while the maximum for the lowest heating rate lingers for some minutes. This effect can be explained by taking into account the differences between the temperature profiles corresponding to different heating rates. As the heating rate increases, the temperature profiles are more pronounced, and the expansion of the lower part of the sample can occur simultaneously with the surface regression due to the pyrolysis. This expansion affects a smaller fraction of the sample before it starts to decompose and the volume to decrease. For very low heating rates, the temperature differences across the sample thickness are lower and the pyrolysis starts after the whole sample has reached the melting temperatures. For heating rate 4 the theoretical thickness at 21 min is lower than 5 mm, confirming this explanation. Figure 10 shows the variation of the sample thickness and the surface and bottom temperatures for heating curve 3. As can be observed, the thickness increases as the temperature increases, and melting occurs over around 30 min. During this time, the sample thickness increases slowly up to a maximum value which coincides with the complete melting of the sample, as can be observed at the end of the temperature plateau of the

442

Energy & Fuels, Vol. 16, No. 2, 2002

bottom surface. Once the thickness starts to decrease, an inflection point in the surface temperature is observed, due to the decrease of the mass of that node. The pyrolysis continues until the complete consumption of the sample. Figure 5 shows the comparison between experimental and theoretical solid conversion for heating rates 2, 3, and 4. Since the maximum value of conversion obtained for heating rate 1 is lower than 0.1, this is not considered. When comparing the temperature profiles, it can be observed that the maximum conversion is reached when the temperature of the rear side of the sample reaches its maximum value. This corroborates the assumption of uniform incident heat flux over the surface, so that no temperature profiles perpendicular to the heat flux exist. When comparing experimental and predicted conversion values, although similar trends are obtained, slight discrepancies do exist, the theoretical conversion being overpredicted for heat flux 3 and 4. Since the temperatures are correctly predicted, an insight into the kinetics of thermal decomposition must be considered. As stated previously, the kinetic equations for this process have been obtained previously by thermogravimetry, using small samples (≈ 3 mg) in which no internal temperature profiles and resistance

Ceamanos et al.

to mass transfer exist.11 The pyrolysis rate obtained under these conditions must be faster than that governing the weight loss in a large sample. Conclusions An experimental and theoretical study of the melting and thermal decomposition of high density polyethylene has been carried out. The mathematical model developed, which includes the melting process and the volume variation adequately predicts the temperature profiles for different heating rates. Some differences are observed between experimental and theoretical solid conversion. The prediction of solid conversion using a slower reaction rate obtained by isothermal thermogravimetric experiments shows much better agreement. This indicates that, given the sample thickness and high viscosity of the melted polymer, an important limitation to mass transfer exists, so the emission of the pyrolysis products is delayed, thus decreasing the solid conversion. Acknowledgment. The authors express their gratitude to C.I.C.Y.T. (Project QUI98-0669) for providing financial support for this work. EF010179E