Hot-pressing process modeling for medium density fiberboard (MDF)

In this paper we present a numerical solution for the mathematical modeling of the hot-pressing process applied to medium density fiberboard. The model is based in the work of Humphrey[82], Humphrey and Bolton[89] and Carvalho and Costa[98], with some modifications and extensions in order to take into account mainly the convective effects on the phase change term and also a conservative numerical treatment of the resulting system of partial differential equations.

Hot pressing is the process in which a mattress composed of wood fibers and resin is cured by applying heat and pressure in a press (see figure 1). Continuum and batch presses do exist, and one of the main issues in reducing the cost of the final product is to reduce the press cycle time. In order to improve the heat transfer between the press platens and the inner layers, some amount of water is added to the mat. Another issue is to adjust the parameters and the temperature history of the cycle in order to obtain a given density profile in the board. Normally, it is desirable to have lower densities at the center of the board in order to increase the mechanical rigidity for a given total mass per unit area. Predicting the influence of these parameters, namely water content, press cycle duration and history (pressure and temperature) is one of the main concerns of numerical models. Many numerical models have been reported to help in predicting the influence of the process parameters in the final product. Among the most complete, we can find the finite difference 2D (axisymmetrical) model presented by Humphrey [1982] and, more recently the 3D model of Carvalho and Costa [1998]. Both consider conduction, phase change of water from the adsorbed to the vapor state and convection. The stress development and the determination of the density profile are not included in these models.
In this paper we present a numerical model which includes all these features and makes some correction to the energy balance equation, as presented in Carvalho and Costa [1998]. The model is based on the finite element model, so that it allows for a more versatile definition of geometry, dimensionality and, eventually, coupling with other packages. It also will allow the use of adaptive refinement, which may be an important issue at the lateral borders, where the hot steam flows from the board to the ambient. However, this issue is not considered in this paper.

Multiphase model
In order to avoid modeling the material down to the scale of the microstructure (the fibers in this case), non homogeneous materials are solved via "averaged equations" so that the intricate microstructure results in a continuum with averaged properties. The averaged equations and properties can be deduced in a rigorous way through the theory of mixtures and averaging operators [Whitaker, 1980].

Energy balance
We will not enter in the details of all the derivations but only for the averaged energy balance equation, which can be found in Appendix B. The referred equation is where T is temperature, t time, C p specific heat, k thermal conductivity, ∇ the gradient operator, ρ s the density of the dry board (solid phase), ρ v vapor density, V g the volume averaged gas velocity, i.e.
where v g is the velocity averaged on the phase, see Appendix B), C pv specific heat of vapor,ṁ evaporation rate, λ latent heat of vaporization of free water, and Q l adsorption heat. The main difference between this equation and that presented by both Carvalho and Costa [1998] and Humphrey [1982] is in the addition of the water evaporation heat term in the convection term instead of considering the phase change effect only on the temporal term. This term should be included because both phases, the solid material and the vapor are in relative motion and we think that its influence is not negligible in a high temperature process. Carvalho and Costa [1998] proposed the following steam mass conservation equatioṅ

Steam mass balance
where MM w is the molecular weight of water, R the gas constant, ǫ board porosity, D v diffusivity of water vapor in the air/vapor mixture and P v vapor partial pressure. Considering that the steam is treated as an ideal gas, then so it may be written, assuming ǫ constant, aṡ This expression is preferable to (3) because it is written in a conservative form that is more agreeable for a numerical treatment. The left hand side term represents the mass interfacial transport and those in the right hand side take into account the mass diffusion and the mass convection. However, it should be noted that this last expression does not have a temporal term as every consistent balance equation does. For example, if evaporation is not considered, then (5) is valid only for a steady situation, which is not in general the case. Then, we rewrite the steam mass balance as This is another difference between our model and that proposed by Carvalho and Costa [1998].

Gas mixture mass balance
Finally, because the gas phase is composed of two main constituents, steam and air, we may use an additional equation for the mass transport of the whole gas phase. Carvalho and Costa [1998] considered where P is the pressure of the gas phase and K g the board permeability tensor. Again, assuming ideal gas law as the state equation for this phase, we arrive to In order to close the system of equations we need to introduce a relationship betweenṁ and (∂P v /∂t). Consider the steam mass balance (6) and the relation that represents the fact that the board moisture content H is composed of vapor and bound water ρ L . If we assume that no liquid phase is considered, then bound water may be transferred to the gas phase only (solid to steam). Soṁ and then The air mass balance equation can be obtained by substracting (6) from (9) Due to the fact that the mean macroscopic diffusive fluxes should be null the air mass balance equation is transformed in the following expression which is very similar to (6) but here valid for the air. Obviously, the air transport equation has no evaporation term.

Summary of equations and boundary conditions
In order to clarify the mathematical model that is finally used for the simulation of hot-pressing process we present the following brief summary of partial differential equations.
• Energy balance equation • Water content balance equation • Air mass balance equation The boundary conditions are the following: • At the press platen: T = T platen (t), air/water mixture in equilibrium with platen temperature V g ·n = 0, no mass flow across the platen ∂ρ a ∂z = 0, no air diffusion across the platen ∂ρ v ∂z = 0, no vapor diffusion across the platen.
• At the center line (r = 0), axial symmetry for all variables • At the mid plane (z = 0), symmetry for all variables • At the exit boundary (r = R ext ), ∂T ∂r = 0, null diffusive heat flux P v = P v,atm , equil. with external air/water mixture, P a = P a,atm , equil. with external air/water mixture.
3 Numerical method The above system of equations contains three main unknowns, the temperature, the moisture content and the air density representing the dependent variables of the problem also called the state variable.
In this work we have used as independent variables the time and two spatial coordinates (3D problems may be computed much in the same way). Due to the physical and geometrical inherent complexity of this problem this may be computed only by numerical methods. For the spatial discretization we have employed finite elements with multilinear elements for all the unknowns. Due to the high convective effects the numerical scheme was stabilized with the SUPG (for "Streamline Upwind -Petrov Galerkin") formulation (see Brooks and Hughes [1982]), otherwise spurious oscillations arise. Once the spatial discretization is performed the partial differential system of equations is transformed into an ordinary differential system of equations likė where U is the state vector containing the three unknowns in each node of the whole mesh. So, the system dimension is 3N where N is the number of nodes in the mesh. The numerical procedure is as follow: Knowing the state vector at the current where j represent an specific node in the mesh. To get the residual, right hand side of (23) the following steps should be done: • Obtain air pressure from the gas state equation ((T, ρ a ) → P a ).
• Obtain relative humidity from sorption isotherms • Obtain vapor density in air from vapor state equation • Compute gradients of T , P from nodal values at the Gauss points using the finite element interpolation. • Assemble the element residual contributions in a global vector.
Once the whole residual vector at t = t n is computed the unknowns variables at the next time step is updated with This kind of scheme, called explicit integration in time is very simple to be implemented but it has two major drawbacks, one is the limitation of the time step to ensure numerical stability and the other is the bad convergence rate for ill-conditioned system of equations. In this application the last disadvantage is very restrictive because the characteristic times of each equation are very different. To circumvent this drawback we have implemented an implicit numerical scheme where the residue is computed at the new state variable U (t n+1 ) instead of using the current value U (t n ). The non-linearities and the time dependency of the state vector makes this implementation more difficult and time consuming but more stable. This non-linear equation in U n+1 is solved by the Newton method, which requires the computation of the Jacobian In order to avoid the explicit computation of the Jacobian, we compute an approximate one by finite differences.
This can be done element-wise, so the cost involved is proportional to the number of elements in the mesh. In our application and despite of the ill-condition of the problem we have found a good convergence at each time step, in average 4 iterations per time step to reduce 10 orders of magnitude in the global residue.
4 Physical and transport properties 4.1 Thermal conductivity Following Humphrey [1982] the influence of board density, moisture content and temperature on the thermal conductivity is considered as an independent correction factor obtained experimentally where κ z is the thermal conductivity in the pressing direction in W/m • K and ρ s is the oven dry density of the material in Kg/m 3 . The moisture correction factor F κ,H (Kollmann and Malmquist [1956] and Humphrey [1982]) assumes H, the moisture content of the board material, in %, and the temperature correction factor F κ,T (Kuhlmann [1962] and Humphrey [1982]) assumes T in • C.

Heat flux direction correction
By far the greater part of conductive heat translation takes place in the vertical plane. However the energy lost from the mattress is largely the result of radial vapor migration from the center toward the atmosphere. The associated horizontal relative humidity gradient lead to a horizontal temperature gradient. Even though this gradient is always lower than the vertical one its influence should be taken into account if multidimensional analysis is required. Ward and Skaar [1963] made experimental measurements and they observed that at a first glance a factor of approximately 1.5 may be a good initial guess before doing some extra measurements. Then

Permeability
The evaporation and condensation of water changes the vapor density and consequently its partial pressure in the voids within the composite.
A vertical pressure gradient leads to the flow of water vapor from the press platens toward the central plane of the board. At the same time an horizontal vapor flow is set up in response to the pressure gradient established in the same direction. The relation between the pressure gradient and the flow features may be assigned to the material permeability. Permeability is a measure of the ease with which a fluid may flow through a porous medium under the influence of a given pressure gradient. Different mechanisms may be involved in this flow, a viscous laminar flow, a turbulent flow and a slip or Knudsen flow.
In this study only the first type is included with the assumption that Darcy's law is obeyed. This may be written as where ∇p is the driving force and v g is the flow variable. K g /µ g is called the superficial permeability with µ g the dynamic viscosity of the gas phase and K g the specific permeability.
The kinetic theory of gases suggests that at normal pressure the viscosity is independent of pressure and it varies as the square root of the absolute temperature, µ ∝ √ T 10 3 < p < 10 6 .
Corrections with the absolute temperature are often considered by Sutherland law with B and C characteristic constants of the gas or vapor with µ in Kg m −1 s −1 . Values for B and C are available from Keenan and Keyes [1966]. For this application the following expression was adopted, µ = 1.112 × 10 −5 × (T + 273.15) 1.5 (T + 3211.0) , with T in • C.

Variation of vertical permeability with board material density
Even though we consider that the press is closed and the density profile is set up and fixed it is included here some conclusions from Humphrey [1982] with results obtained by Denisov et al. [1975] for 19 mm boards and others from Sokunbi [1978]. The data to be fitted are the following

Horizontal permeability
Sokunbi measures included in figure 2.7 of Humphrey [1982] shows the relation between the board thickness in mm with horizontal permeability. For approximately 15mm board thickness and ρ s = 586Kg/m 3 the horizontal permeability is 59 times the vertical value, in agreement with the values assumed by Carvalho and Costa [1998].

Steam in air diffusivity
The interdiffusion coefficient of steam in air can be calculated from the following semi-empirical equation [Stanish et al., 1986], where the diffusivity is in m 2 /sec, pressure in N/m 2 and T in • K.

Vapor Density
For the pressure range likely to occur during hot pressing (between 10 3 and 3 × 10 5 N/m 2 ) a linear relationship between saturated vapor pressure and vapor density may be assumed. Fitting experimental data Humphrey [1982] proposed the following expression with ρ v in Kg/m 3 , P sat in N/m 2 and the relative humidity HR in %. This can be deduced from the relative humidity definition and applying ideal gas law for gaseous phase TakingR = 8314J/Kmol/K and MM w = 18Kg/Kmol with T ≈ 360 • K we obtain (35).

Saturated vapor pressure
Following the Kirchoff expression with data presented by Keenan and Keyes [1966] we include here the following equation, log 10 P sat = 10.745 − (2141.0/(T + 273.15)), with P sat in N/m 2 and T in • C.

Latent heat of evaporation and heat of wetting
Using Clausius-Clapeyron equation in differential form and after some simplification the latent heat of vaporization of free water may be written as λ = 2.511 × 10 6 − 2.48 × 10 3 T, with λ in J/Kg and T in • C. For the differential heat of sorption we follow Humphrey [1982] that used (see Bramhall [1979]) with Q l in J/Kg and H in %.

Specific heat of mattress material
It is computed by adding the specific heat of dry wood and that of water according to the material porosity and assuming full saturation. The specific heat of dry mattress material is taken as 1357J/Kg/K and the specific heat of water has been taken to be 4190J/Kg/K. From Siau [1984] the expression for specific heat of moist wood is T in • K and H in %.

Porosity
According to Humphrey [1982] the volume of voids within the region may be computed with where ǫ is the porosity, ρ is the density of the region and ρ s is the dry density of the board material.
In Carvalho and Costa [1998] they included the expression from Suzuki and Kato [1989] ǫ = 1 − ρ s 1/ρ f + y r /ρ r 1 + y r , where ρ r is the cure resin density, ρ f is the oven dry fiber density and y r is the resin content (resin weight/board weight). In Carvalho and Costa [1998], the authors used y r = 8.5% with ρ f = 900Kg/m 3 and ρ r = 1100Kg/m 3 .

Numerical results
In this section we present some results that can be compared with those reported in Humphrey [1982]. This numerical experiment allows the validation of the mathematical model and its numerical implementation for future applications to hot-pressing process simulation and control. This experiment consists of a round fiberboard of 15 mm of thickness and 0.2828 m of radius that according to its axisymmetrical geometry needs as spatial coordinates only the radius r and the axial coordinate z assuming no variation in the circumferential coordinate. The axisymmetrical domain is discretized in 20 × 20 elements in each direction with a grading toward the press platen and the external radius as may be visualized in figure 2. In order to follow the same assumptions as in that work we fixed the air density to a very low value (ρ a ≈ 10 −6 ), as if the press was close with no air inside the fiberboard. The boundary conditions are as in section 2. For the press platen temperature we have applied a ramp from 30 • C at t = 0 to 160 • C at t = 72 seconds with a least square fitting from data in Humphrey [1982]. As initial conditions we have assumed a uniform temperature of T (t = 0) = 30 • C in the whole solid material with a uniform moisture content of H = 11%. The external atmosphere was considered to be at T atm = 30 • C, HR atm = 65%, in such a way that the internal moisture content is, at the initial state, in equilibrium with the external atmosphere.
In the next section we include the results obtained by using the model above cited to the original problem presented in Humphrey [1982]. Next we present some other experiments showing some phenomena that deserve more attention for futue studies. Figure 3 shows the temperature distribution in time and with the axial coordinate at r = 0 (centerline). We can note the penetration in the axial direction of the temperature profile in time for t = 1, 10, 50, 100, 200, 300 and 400 seconds. Figure 4 shows the same kind of plot for moisture content. For r not to close to the external radius, the problem is almost one dimensional in the z direction. As the thermal front penetrates into the board water evaporates. This vapor advances to lower pressure regions near the symmetry plane and, as it encounters lower temperatures, it condenses releasing heat. This process can be clearly seen from the wave in moisture content (figure 4) exceeding the initial water content of 11% and results in an improvement in the heat transfer with respect to the pure conduction case. Also this phenomenon is responsible of the change in curvature of thee temperature curves, mainly at t = 10, 50, and 100 sec (see figure 3). The total water content in the board at a given instant can be found by integrating the bound water content and the water in vapor phase. However, this last is negligible. We can see in figure 4 that the depression in water content near the board (for instance at t = 400sec), is larger that the water enrichment in near the center plane. This is due to water migration from the center of the board to the external radius, where it flows to the external atmosphere. The following figures show similar plots but at different locations,  Figure 11 shows a three-dimensional view for moisture content represented by the third coordinate axis at t = 200 seconds as a function of r, z.

Original numerical experiment
These results are in good agreement with those presented by Humphrey [1982]. However the cited author did not present his results at some locations that in our opinion should be treated with some care, for example at the external radius.

Further numerical experiments
In Humphrey [1982], results for moisture and temperature in the vertical and radial directions at both central planes, r = 0 and z = 0 respectively are included. No mention about the vertical distribution at r = R = 0.2828m or about the moisture content at press platen. Our results present some overshooting in the temperature profile very close to the radial exit contour and at the first moment we though about a spurious numerical problem, but it is due to large variations of the magnitude of vapor pressure and density at the boundary. In typical runs, vapor pressure varies from near 2 atm in the center of the board to 0.01 atm at the external radius. We think that this problem will be fixed if we solve for the air density also, but then a very fine grid will be necessary at the exit boundary, since large variations of the vapor molar fraction are expected. (see figure 12). Molar fraction varies from nearly 1 at the interior of the mat to a 2% at the external atmosphere. This variation is produced in a thin layer of width δ proportional to the diffusivity of vapor in air which is very small.

Conclusion
We presented a numerical model for the heat and mass transfer in the hot-pressing model of a MDF fiberboard. The model includes convective effects on the phase change term and also a conservative numerical treatment of the resulting system of partial differential equations. Convective effects are responsible of an increase in heat transfer from the mat Total pressure Vapour partial pressure air δ Figure 12: Boundary layer in vapor partial pressure at external radius platen to the center of the board due to water vapor evaporation and condensation. Two-dimensional simulations allow to estimate border effects.
where h is enthalpy. For the other phases (solid and bound water) a similar expression holds, but neglecting the advective term. Applying the volume average operator [Whitaker, 1980], we arrive to the following equation averaged on the gas phase ǫ g is the volumetric fraction of phase g (i.e. gas) and X g is the average of quantity X on the volume occupied by phase g The term Q g is the total enthalpy flux through the solid-gas interface Γ where (X) g is the value of property X on the g side of the interface and w is the velocity of the interface. Assuming that h g is constant on all Ω g (for a certain volume control) then whereṁ is the rate of mass of water being evaporated. A common drawback of averaged equations is that, when products of variables like ρh appear in the microscopic equation, the average of the product ρh g is obtained in the averaged equation. Now, it is not true that so that the averaged equation contains more unknowns than the original equation. A common assumption is that no correlation exists between variables and so (49) is approximmately valid. This can be justified, for instance, if the variations of each quantity around the mean is small. Then, applying the volume average operator over the gas, solid, and bound water phases and assuming no correlations between variables we obtain the following averaged equations for the three phases ∂ ∂t (ǫ g ρ g h g ) + ∇ · (ǫ g ρ g v g h g ) = ∇ · (ǫ g k∇T g ) −ṁh g gas phase, ∂ ∂t (ǫ s ρ s h s ) = ∇ · (ǫ s k∇T s ) solid phase and , ∂ ∂t (ǫ l ρ l h l ) =ṁ h l liquid phase.
Also, the average operator is dropped from here on, and a subindex g or s implies averaging on that phase. Also, V g the volume averaged gas velocity is Note that in the body of the text V g is used instead of vg. Now, summation of these three equations gives ∂ ∂t (ǫ g ρ g h g + ǫ s ρ s h s + ǫ l ρ l h l ) + ∇ · (ǫ g v g ρ g h g ) = ∇(ǫ g k∇T g + ǫ s k∇T s ) +ṁ (h l − h g ). (52) Now, where k eff is the average conductivity of the solid+water+gas mixture. The gas is assumed as an ideal mixture, so that the enthalpy is the sum of the enthalpy of its constituents, and neglect the contribution of the air constituent so that Taking a reference state for the entalphy at a point on the adsorbed state We also neglect the entalphy of the gas phase with respect to the solid+water phases and put where ρ eff , and C p eff are averaged properties for the moist board, as a function of temperature and moisture content. Finally, the averaged equation is This is equivalent to (1) through relation (2) and assuming T ref = 0 • C.