Collection Mémoires et thèses électroniques
Accueil À propos Nous joindre

Chapitre 5 Modeling of heat and mass transfer during MDF hot pressing: Part 2 Finite element modeling of batch and continuous pressing

Table des matières

Un modèle mathématique pour le pressage à chaud des panneaux MDF qui peut être appliqué aux procédés en lot et en continu a été dans une publication précédente. On présente une stratégie de solution numérique basée sur ce modèle. On emploie la méthode des éléments finis pour la simulation bi et tridimensionnelle de ce modèle. Le système fortement non linéaire des équations gouvernantes du procédé de pressage à chaud est résolu en utilisant un programme local nommé MEF++. Le modèle prédit l’évolution dans le temps des variables telles que la pression gazeuse, la teneur en humidité et la température durant le procédé. Ce système de trois équations couplées représente une description physique complète du pressage à chaud pour un panneau MDF fabriqué au niveau macroscopique. On utilise le schéma d’Euler implicite pour la discrétisation en temps et la méthode de Newton pour linéariser le système. Les résultats numériques présentés sont en accord avec les résultats expérimentaux obtenus au laboratoire et en industrie. Cela démontre le potentiel de la méthode des éléments finis utilisée comme outil d’analyse performante pour la fabrication des panneaux MDF dans un contexte industriel.

Mots clés : pressage à chaud des panneaux MDF, méthodes des éléments finis, transfert de masse et de chaleur, équations couplées

A mathematical model of heat and mass transfer in the mat for MDF hot-pressing applicable to batch and continuous processes was previously introduced. A numerical solution strategy based on this model is now presented. The numerical technique is based on a finite element method and can be used for both two and three dimensional simulations. The highly non-linear system of equations governing the hot-pressing process is solved using an in-house program called MEF++. The model predicts the time evolution of variables such as total gas pressure, moisture content and temperature during the complete process. This system of three-coupled equations represents a complete physical description of the hot-pressing for randomly formed MDF panel at the macroscopic level. The solution procedure uses a fully implicit time-stepping scheme and Newton’s method to take care of the non-linearities. Numerical results are in relatively good agreement with laboratory and industrial results. This shows the potential of the finite element method as a tool for the design of MDF panels in an industrial context.

Key words: MDF hot-pressing, finite element method, mass and heat transfer, coupled equations

Modeling the hot pressing of wood-based composite panels is of great interest for industry. It has been the object of several research efforts, where various approaches were used. Cylindrical board geometry was considered by Humphrey and Bolton (1989) for the modeling of medium density fiberboard (MDF) hot-pressing. This kind of geometry is still used by some researchers to simplify the boundary conditions (Carvalho and Costa 2003). In realistic conditions, rectangular panels have to be considered (Suo and Bowyer 1994, Zombori 2001, Thoemen 2000, García 2002) and simulations must be performed in two or three dimensions. Only Carvalho and Costa (1997) present comparisons between two and three dimensional models for heat and mass transfer during MDF hot-pressing. Clearly, two-dimensional calculations are less expansive but this simplification is not always possible. It is thus important to develop a numerical tool allowing for three dimensional simulations in various geometries.

Different numerical methods were used to solve models of heat and mass transfer applied to wood-based panels hot pressing including the finite element method (Hubert and Dai 1999, Nigro and Storti 2001), the finite difference method (Harless et al. 1987, Suo and Bowyer 1994, Humphrey and Bolton 1989, Thöemen 2000, Zombori 2001, Carvalho 2003) and the control volume method (Turner and Perré 1995, García 2002). The choice of the most appropriate numerical solution strategy required to obtain good results in the shortest computing time possible is not an easy task. Many numerical simulations have been performed for wood-drying problems such as the work done by Turner and Perré (1995). There is however a lack of information concerning the numerical simulation of the hot-pressing process of wood-based panels and there is a definitive industrial need for an efficient software in this field.

Taking into account previous work done on predictive models and knowing the advantages of the finite element method, the main objective of this study was to implement a new physical model of heat and mass transfer in the MDF mat during hot-pressing problems and validate it by comparing the simulations with laboratory and industrial results. The finite element method was chosen since it allows a more flexible definition of the geometry and a more accurate definition of complex boundary conditions.

Neglecting the contribution of the source terms in equation (5-3), we multiply the first equation by a test function , the second by and, the third by . We integrate over the domain, using the divergence theorem (Reddy, 1993) to get:

where is the exterior unit vector of the boundary.

Boundary conditions must also be provided. From the above variational formulation, it is possible to impose , , and on a given section of the boundary . This is called a Dirichlet or essential boundary condition. On the other hand, natural or Neumann boundary conditions, arising from the integration by parts can also be applied . Only one of the two boundary conditions (essential or natural) can be specified. From the surface integral terms of Equations (5-9), (5-10) and (5-11), one has to impose on the boundary:

For hot-pressing problems, it is common to impose , and on the sections of of the boundary (and consequently on ). Most importantly, “hot-pressing conditions” expressing heat transfer from the hot-platen to the mat have to be imposed on the section of the boundary:


: Mass convection coefficient based on pressure (P) in air [kg m-2s-1Pa-1]

: Mass convection coefficient based on moisture content (M) in wood [kg m-2s-1%-1]

: Thermal convection coefficient [Wm-2s-1K-1]

: Some external gas pressure [Pa]

: Some external moisture content [fractional]

: Some external temperature [K]

Using (5-15), (5-16) and (5-17), the variational formulation becomes:

It is important to note that:

  • This is a model with three unknown state variables (P, M, and T). Every other variable must be expressed in terms of these three variables.

  • The problem is highly non linear.

  • Many material properties have to be provided and in several instances, little is known about them for wood based panels in general and MDF in particular.

A fully implicit Euler scheme is employed to discretize the time derivatives. Starting from the solution at time ( , and being given as initial condition) this gives:

At each time step, a non linear system is obtained and solved by Newton’s method. Linear and quadratic finite elements can be used in the program. Quadratic elements were used in two dimensions while linear elements were used in three dimensions to reduce the computational burden.

Two case studies were considered: laboratory batch pressing and industrial continuous pressing. These two cases will be presented below.

To compare numerical results with experimental data, we used a homogeneous MDF mat made from 90 % black spruce ( Picea mariana (Mill.) BPS ) and 10 % balsam fir ( Abies balsamea (L.) Mill ) fibers obtained from Uniboard Canada Inc., Panneaux MDF La-Baie plant . The initial mat moisture content was 9.3 %. The press temperature was constant at 185 ºC and the target board density was 710 kg m-3. The final thickness of the mat was 11 mm and it was manufactured using a 600 × 600 mm2 laboratory press. More details about the MDF panels manufacturing parameters are given in Table 5-1.

Mat temperature was measured using K-type thermocouples positioned at the surface and center of the mat respectively. Gas pressure was measured at the mat center using a PressMAN gas pressure transducer (Figure 5-1). The press closing time was between 80 and 90 s at a maximum pressure of about 7 MPa to achieve the final thickness. The pressure was then reduced to 5.5 MPa and kept constant for 177 s. Finally, the press hydraulic pressure was gradually reduced to zero and the press opened within approximately 60 s for a total pressing cycle of approximately 330 s. Position control of mat thickness was determined by a LVDT on the top platen.

Considering the same plane geometry described in Figure 5-2, the boundary conditions specified on the 3-D industrial mat sections are a little more complicated. On top and bottom surfaces, we have to distinguish the parts of the mat inside and outside the press. If is the feed speed [m s-1] as indicated in Figure 4 and t is the pressing cycle time [s], then at time t , the press position is defined as:

It is important to note that after 50 s the mat is completely inside the press. The part of the mat for which is still outside the press. On this part, we suppose that we have convective boundary conditions for both the gas pressure and the moisture content. Inside the press, we impose vanishing Neumann conditions. More precisely, the boundary conditions are the following:

On top surface and bottom surfaces:


Boundary conditions on edges 1 and 2:

Boundary conditions on edges 3 and 4:

where [K] and [Pa].

Other boundary conditions can be easily imposed but the above set was perceived as the most representative to the industrial process.

Only one mesh was used in that case consisting of a 16x16x16 = 4096 linear elements. The details of the input data for laboratory and industrial MDF panels pressing are listed in Table 5-3. The time step was set to 2 s. The numerical results show the evolution of temperature, gas pressure, and moisture content during hot-pressing.

Numerical results obtained for moisture content as a function of mesh size and dimensions considered (2-D or 3-D) are shown in Figure 5-5. It is obvious that two and three dimensional simulations give essentially the same results. This is no surprise for such a simple geometry. Mesh refinement has no significant impact on the results for surface moisture content. We can observe a small reduction of surface moisture content from 2-D to 3-D models but this difference can also be attributed to the mesh size.

On the other hand, the evolution of moisture content is in agreement with our expectations. Surface moisture content decreases more quickly than center moisture content. This is due to vapor flow towards the interior of the mat.

The temperatures at the surface and center calculated by the numerical model follow a pattern similar to the measured data (Figure 5-6). The surface temperature reaches 160 ºC quickly due to the direct contact with the platen and then gradually rises to the constant platen temperature as expected. However we note that the temperature at the center starts to rise sooner in the numerical simulations. A faster heat transfer to the center by heat conduction or heat convection by gas bulk flow calculated by the model can explain this discrepancy. This can be due to improper mat thermal conductivity or gas permeability used in the model. Also, since the simulation is performed using the final mat thickness, the center temperature starts to rise sooner than in the experimental case due to the shorter distance involved in the heat transfer. Calculated and measured temperatures are close to press temperature at the end of hot pressing at 170 ºC and 180 ºC (300 s) respectively. Temperatures predicted at surface and center of the mat are in relatively good agreement with experimental data (Figure 5-6). Center temperature reached 154 ºC at 200 s while it was about 146 ºC for the measured data at the same time.

Figure 5-7 shows a typical temperature profile developed in the thickness direction. It has a typical U shape: a high temperature occurs in the external layers very close to the temperature of the press and a lower temperature in the core layer (103 °C).

As we did not have experimental data on the moisture content of MDF panels for these specific conditions, comparisons can be made with the Malmquist (1959) wood sorption model to predict moisture content through mat thickness (Figure 5-8). The model was used to calculate equilibrium moisture content corresponding to temperature and equilibrium moisture content in the mat. Details on this mathematical model to predict the equilibrium moisture content can be found in Vidal and Cloutier (2005). It can be seen that the numerical results at 300 s are lower for wood composites (more or less 2 %) than for wood (more or less 6 %). It is well known that the equilibrium moisture content for wood composites is lower than for solid wood (Wood Handbook, FPL 1999). Considering this difference, numerical results for moisture content of MDF panels present similar trends during the hot-pressing such as results presented by García (2002). We can therefore consider the numerical results as acceptable. The effect of different oven-dry density is also presented in Figure 5-8. These results are also in agreement with García (2002): with denser wood it is more difficult to extract the bound water from wood cell walls. Higher density panels lead to higher moisture content values.

Experimental and numerical results for gas pressure are displayed in Figure 5-9. Since we know the mat thickness variation from experimental measurements (also in Figure 5-9), we can distinguish three stages: press closing, mat consolidation and polymerization. This variation of mat thickness during consolidation was not considered in the current numerical model as for other models proposed in the literature. This is clearly a limitation of these models.

As expected, gas pressure increases gradually above atmospheric pressure. However, the trends in the numerical prediction are significantly different from what is observed in the experimental data. According to Figure 5-9, no change is detected at the center of the mat by the PressMAN probe for about 70 s. At this time, a small pressure increase is observed at the end of mat consolidation. On the other hand, the numerical results show a rapid pressure increase until a plateau is reached. The pressure value corresponding to the plateau is about 135 kPa, significantly smaller than the final value in the experimental case. The differences between numerical and experimental results can probably be attributed to different sources but the most important is undoubtedly the absence of thickness variation in the numerical model. Other sources of error are the sorption model chosen, inadequate boundary conditions and inaccurate values of the convective coefficients.

A 3-D simulation was performed for the case of industrial continuous pressing of MDF panels using a mesh of 4096 elements. The experimental data obtained from the measurements described in Figure 5-4 are displayed in Figures 5-10 and 5-11 and compared with the numerical results. Figure 5-10 shows measured and calculated center temperature. The calculated temperature is in relatively good agreement with the measured values. Although the pattern of the curves is slightly different, the calculated center temperature reached 99 ºC, very close to the measured value of 111 ºC at 364 s. The calculated center temperature does not show the plateau observed for experimental data in the first 10 m of the press.

The calculated gas pressure is in good agreement with the experimental data as can be seen in Figure 5-11. The agreement with experimental data is much better in this case compared to first case. A possible explanation is the weaker influence of edge boundary conditions in the second case due to the dimensions of the mat considered. It must be emphasized that no adjustment factor was used in numerical results for both cases analyzed.

Gas pressure at point 2 (Figure 5-4) at 13 % of the width is consistent with the results of Thoemen (2000) where gas pressure was measured by inserting probes at 40 % and 22 % of the mat width (Figure 5-12). The general trend is similar in the three cases. As we approach the mat edges, the gas pressure decreases to reach equilibrium with the atmospheric pressure because vapor escapes continuously by the edges of the mat.

The gas pressure evolution for the second case is presented in Figure 5-13. At t = 2 s, the mat enters into the press and the compression starts. The gas pressure is higher in the first zone of the mat at this time. The remaining of the mat is still out of the press and thus at lower gas pressure (1.0132*105 Pa) since it takes 50 s for the 5-m long mat considered to completely enter into the press. At t = 100 s, the mat presents a gas pressure distribution in three zones: higher pressure for the first zone (1.1895*105 Pa), transition pressure for the center zone (1.1763*105 Pa) and lower pressure for the third zone (1.1598*105 Pa). At t = 300 s, the gas pressure is concentrated in the center and the third zone of the mat with values of 1.2806*105Pa and 1.2841*105 Pa respectively. This pressure is extended to the length of the mat. At 370 s the higher gas pressure is still concentrated in the same zones of the mat with slightly higher values (1.2989*105 and 1.3049*105 Pa). From Figure 5-13, we can observe lower values of pressure at the four edges of the panel as opposed to the center. During hot-pressing, gases are transported in the thickness direction but once inside the mat, they are forced towards the open edges thus developing a pressure gradient from the center to the edges of the panel. When the temperature increases, an important generation of vapor is produced. This phenomenon causes a pressure increase especially at the center of the mat. This effect was observed during hot pressing because the vapor is continuously forced out by the edges of the mat. This is confirmed by the numerical results obtained from Figure 5-13.

In Figure 5-14, surface moisture content (point 1 from Figure 5-4) presents a similar behavior as for the first case. The hot-platen produces a quick drop of moisture content at the surface of the mat. Moreover as can be seen in Figure 5-15, a non-homogeneous distribution of moisture content is observed on the top surface of the mat. At t = 2 s, compression and heating of the mat start and the loss of moisture is immediate in zone 1 (2%, 4%, 6% and 7.9%). The remaining of the mat is still at the initial moisture content. At t = 100 s, the moisture content is lower in the first region (1.3 %) and increases slightly towards the end zone (2.9 %). A strong moisture content gradient starts to develop at the edges of the mat due to the boundary conditions. In industrial hot pressing, moisture content is allowed to escape through the lateral edges only. This underlines the influence of the convective parameter h M_2 which was most likely too small. This value is unfortunately not well known. A larger value would probably produce a smaller moisture content gradient. At t = 300 s, the moisture content shows a linear progression along mat length (zone 1: 1.32%, zone 2: 1.79 % and zone 3: 2.4 %) but on the edges, the boundary layer gets stronger. Finally, at t = 370 s, the general pattern is similar. In the four corners, oscillations are observed due to the incompatibility of the boundary conditions.

The objective of this study was to demonstrate the potential of the finite element method to predict the behavior of the most important variables affecting heat and mass transfer involved in the MDF hot-pressing process in 2-D and 3-D. Despite some discrepancies observed for gas pressure in the laboratory batch pressing case, numerical results are in reasonable agreement with measured data for temperature and moisture content. However, in the industrial continuous pressing case, the predictions are satisfactory for the three state variables. The numerical results show that the model predicts different behaviours for the pressure, temperature and moisture content for the two pressing strategies.

The main conclusion of this study is the lack of accurate data on physical properties of wood composites such as the diffusion D or the convection coefficients hP and hM. These coefficients have a major impact on the numerical results, particularly those appearing in the boundary conditions. In the presented numerical results, no adjusting factor was used to improve the match between measured and calculated results. We have used the most reasonable values found in the literature.

In future work, inverse engineering techniques could be used to improve the parameters values. On the other hand, this model will be eventually coupled with a mechanical model that would take into account the variation of mat thickness during the process.

The authors are grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) for supporting this project under grant no CG 029167. They would like to thank to Uniboard Canada Inc., Panneaux MDF La Baie panel plant for providing the wood fibers and to Pierre Martin for his cooperation in this research. They wish to thank to Serge Plamondon for technical assistance at Laval University.

© Marcia Vidal Bastias, 2006