Collection Mémoires et thèses électroniques | |
Accueil À propos Nous joindre |
Table des matières
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.
For modeling purposes, the mat is considered as a homogeneous material, since morphological characteristics of the fibers were not considered. The model provides a complete description of moisture content, temperature, and total gas pressure profiles in the mat and their evolution with time. The characterization of the heat and mass transfer includes input data such as thermal conductivity and gas permeability. This last property was measured in-situ to represent the characteristics of the MDF panel employed during the numerical simulation. The remaining physical and transport properties have been taken from available literature and adapted as well as possible to the characteristics of the mat. A detailed description of the model of heat and mass transfer in the MDF mat during hot pressing is given in Vidal et al . (2006).
The basic assumptions of our physical model of heat and mass transfer in the MDF mat during hot-pressing are the following:
The wood fibers are assumed to be below the fiber saturation point (FSP). Water in wood was thus assumed to be in the state of water vapor or bound water.
The diffusion mechanism occurring simultaneously with the bulk flow was modeled according to Fick’s first law.
Gas and vapor flows were expressed according to Darcy’s law.
Potential and kinetic energies were neglected in the energy equation.
The total internal-energy density was assumed to be the sum of the internal-energies of the wood (from dry-wood and bound-water), air and vapor components.
The total energy flow in the system is due to heat conduction across the hot-platen interface (the most important source of heat transfer), heat convection in the mat due to gas and vapor bulk flow, and heat loss by convection through mat edges exposed to ambient temperature. Radiation heat transfer is neglected.
Laminar viscous flow is assumed.
Variation of the average mat density was considered in explicit form in the model. It was calculated from the definition of density and press platen position. Mat volume was defined as a function of mat thickness from experimental data obtained from a displacement transducer installed on the press platens.
Hot-pressing involves simultaneous heat and mass transfer into the mat. It is thus necessary to analyze vapor and gas transport through this medium. To achieve this goal, a system of three coupled equations describing the hot-pressing process was used (Vidal et al. 2006). The resulting set of governing equations is:
where:
: Density of gas [kg m^{-3}]
: Density of vapor [kg m^{-3}]
: Density of internal energies [J m^{-3}]
: Time [s]
: Total pressure [Pa]
: Temperature [K]
: Moisture content [fractional]
: Divergence of gas flux vector [kg m^{-3}s^{-1}]
: Divergence of water vapor flux vector [kg m^{-3}s^{-1}]
: Divergence of convection and conduction flux vector
[W m^{-3}]
: Source term for gas [kg m^{-3}s^{-1}]
: Source term for vapor [kg m^{-3}s^{-1}]
: Source term for internal energies [W m^{-3}]
The addition of the diffusion term is one of the differences with the model proposed by García (2002).In its simplest form, MDF hot-pressing can be represented by macroscopic partial differential equations as expressed in equations (5-1), (5-2) and (5-3). Moreover, it was shown in Vidal et al. (2006) that:
where:
We end up with the following system:
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:
where:
: Mass convection coefficient based on pressure (P) in air [kg m^{-2}s^{-1}Pa^{-1}]
: Mass convection coefficient based on moisture content (M) in wood [kg m^{-2}s^{-1}%^{-1}]
: Thermal convection coefficient [Wm^{-2}s^{-1}K^{-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 mm^{2} 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.
Table 5-1 MDF panel manufacturing parameters
Parameter |
Value used |
Board size Mat moisture content Target panel density Wax content Resin type Resin solid content Resin viscosity Resin gel time Resin content Catalyst used pH after catalyst Press platen temperature Total press closing time |
560 × 460 × 11 mm^{3} 9.3 % (based on oven-dry wood) 710 kg m^{-3} 1 % (based on oven-dry wood) UF Commercial 65 % 0.335 Pa s > 15 minutes 14 % (based on oven-dry wood) Solution of NH_{4}CL at 30% 7.0 185 ºC 330 seconds |
The geometry used for the numerical simulation of laboratory hot pressing is presented in Figure 5-2.
The boundary conditions (according to Figure 5-2) are the following:
On top surface and bottom surfaces:
On edges 1, 2, 3 and 4:
where:
[K] and [Pa]
Comparisons with the numerical solutions in 2-D and 3-D will be presented. Mat consolidation and venting were not considered. A simulation run time of 300 s was used and considered as reasonable compared to industrial continuous pressing of MDF panels for a target thickness of 11 mm. Information about the meshes employed for the simulations in two and three dimensions is given in Table 5-2. It is important to note that all of them are homogeneous meshes.
Table 5-2 Characteristics of meshes employed for the first case: laboratory batch pressing.
Figure 5-3 Platen temperature and applied pressure used for the second case: industrial continuous pressing.
A representation of the pressing cycle and platen temperature during continuous pressing is given in Figure 5-3. The continuous press considered is 3 m wide and 37 m long. As can be easily seen, the applied pressure is higher at the beginning of the hot-pressing in order to reach the desired thickness and density into the mat. In addition, platen temperature is not constant and varies with position in the press.
A mat moving into the continuous press is shown in Figure 5-4. The measurements of temperature and gas pressure were made with thermocouples and PressMAN pressure transducers respectively in the center of the mat. They were inserted at 40 cm from the edge (see point 2 in Figure 5-4, i.e. at 13 % of mat width). Point 1 will be used only to report numerical simulation results (Figure 5-14).
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:
where:
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.
Figure 5-5 Numerical results for surface and center moisture content as a function of mesh size and type.
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.
Figure 5-6 Measured and numerical results for temperatures during hot-pressing for first case: laboratory batch pressing.
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-7Calculated two-dimensional temperature profile after 100 s of hot-pressing for first case: laboratory batch pressing.
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).
Figure 5-8 Evolution of moisture content for wood using Malmquist (1959) equation and modeled center moisture content during MDF hot-pressing at different oven-dry densitiesfor first case: laboratory batch pressing.
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.
Figure 5-9 Measured and calculated gas pressure, and mat thickness for first case: laboratory batch pressing.
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.
Figure 5-10 Measured and calculated MDF center temperature for the second case: industrial continuous pressing.
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.
Figure 5-11 Measured and calculated gas pressure for the second case: industrial continuous pressing.
Figure 5-12 Measured gas pressure for different positions into MDF mats. (Results of Thoemen 2000 and current study at 13 % of the width).
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*10^{5} 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*10^{5 }Pa), transition pressure for the center zone (1.1763*10^{5} Pa) and lower pressure for the third zone (1.1598*10^{5} 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*10^{5}Pa and 1.2841*10^{5} 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*10^{5} and 1.3049*10^{5} 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.
Figure 5-13 Distribution of gas pressure (P*10^{5}) in the XZ plane for the second case: industrial continuous pressing.
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.
Figure 5-14 Calculated surface moisture content for the second case: industrial continuous pressing.
Figure 5-15 Distribution of moisture content in the XZ plane for the second case: industrial continuous pressing.
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 h_{P} and h_{M}. 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.
: Thermal coefficient at top/bottom in the panel
: Thermal coefficient at the left/right edge in the panel
: Mass coefficient at top/bottom in the panel for moisture content
: Mass coefficient at left/right in the panel for moisture content
: Mass coefficient at top/bottom in the panel for pressure
: Mass coefficient at left/right in the panel for pressure
© Marcia Vidal Bastias, 2006