samedi 16 août 2008

Numerical analysis on RIM of polyurethane foam


Introduction


Rigid polyurethane foam is widely used at various temperatures as a thermal insulator, a lightweight structure, and a shock absorber in industry. Foam reaction injection molding (FRIM) is one of the most popular and useful processes for producing polyurethane foam whose thermal conductivity, impact strength, elastic modulus, and foam density are the important material parameters. Therefore, it is essential to control the properties during processing in order to obtain the required products. FRIM is a process that consists of many phenomena including mixing, chemical reaction, bubble nucleation, bubble growth, and heat, mass, and momentum transports. Each of these is a complex subject that has been studied intensively and it is a difficult task to consider all of them simultaneously. For this reason, there have been few studies reported on FRIM. Baser and Khakhar modeled the dynamics of a physical blowing agent (CFC-11) and a chemical blowing agent. Water is used as the chemical blowing agent because it reacts with isocyanates to yield carbon dioxide gas. CFC-11 is a chlorofluorocarbon (trichlorofluoromethane) and used as a physical blowing agent. It was assumed that the evaporation of CFC-11 is controlled by heat generation and the carbon dioxide generation is controlled by the water–isocyanate reaction rate. But they did not consider the resin flow or heat loss. Arai et al. carried out experiments in which premixed foam reactants were poured into an L-shaped mold and allowed to fill the mold as the mixture was foamed and expanded. The distribution of foam density was then evaluated in the solid part. They showed that foam density was affected by the pressure variation, although the pressure is not the main factor that controls foam density. Lefebvre and Keunings studied the continuous flow of chemically reactive polymeric liquids in two dimensional geometries using a finite element method. The gelling reaction of polyurethane and the blowing reaction of the chemical blowing agent were taken into account. They claimed that the foam density was a function of temperature, which varied due to the exothermic reaction. Mitani and Hamada predicted flow patterns in the foaming process of polyurethane by considering foam expansion caused by the bubble growth. Three-dimensional control volume finite element method (CVFEM) was used to solve the Stokes equations under isothermal conditions. Density change was measured from the self-expansion of foam by assuming that the density is a function of time only. The numerical prediction of the flow front advancement was compared with the experimental observation. Seo et al. investigated the behavior of a self-expanding fluid numerically with the assumption that the density decreased exponentially with time. The self-expanding fluid showed significantly different flow behaviors from that of the Newtonian fluid of constant density. In this study, foam reaction injection molding was studied numerically by considering chemical reactions, foaming, and mold filling. Theoretical models were developed based on the model proposed by Baser and Khakhar in order to predict the temperature, density, and viscosity variations with respect to time. Ideal mixing and rapid bubble nucleation were assumed. Both gelling and blowing reactions as well as evaporation of physical blowing agents were considered simultaneously. Based on the theoretical model, fully three dimensional numerical simulation of the polyurethane foam flow was carried out to predict velocity and temperature fields, flow front advancements, pressure and density distributions, and thermal conductivity variations. A pressure based finite volume method was selected as the numerical method because fully three-dimensional flow of polyurethane foam had to be simulated by considering gelling and blowing reactions, momentum balance, and heat transfer simultaneously. When the polyurethane foam is used as the insulation material for refrigerators, density and thermal conductivity variations are critical parameters for the performance of the product. Therefore, it is extremely important to understand flow behavior, bubble size distribution, density variation, and thermal conductivity distribution



2. Theoretical model

There are two main processes in the polyurethane foam molding: polymerization and expansion. The polymerization process increases both viscosity and temperature of the foam because the reaction between isocyanate and polyol is highly exothermic. The polyurethane foam is expanded by supplying gas into the nucleated bubbles and ends up with final cellular structure. The gas is usually supplied by physical or chemical blowing agents. In order to model the dynamics of polyurethane foam formation, certain assumptions are needed. Kinetic parameters for the gelling and blowing reactions are assumed to be independent of each other. It is also assumed that the polyurethane foam is a homogeneous phase, the evaporation of the physical blowing agent is controlled by the heat generation, and generation of the carbon dioxide is controlled by the water–isocyanate reaction.

Energy balance

If there is no fluid flow and the material parameters are constant in the temperature range of interest, the energy balance in a control volume under adiabatic condition is expressed as

C is the heat capacity and r is the mass of eachcomponent per unit mass of polymerizing mixture (mixture of un-reacted polyol and isocyanate). The subscript p means the polymerizing mixture, CO2 carbon dioxide, W water, BG physical blowing agent in the gas phase, and BL physical blowing agent in the liquid phase. X is the fractional conversion and (KDH) denotes the heat of reaction, whose subscripts OH and W mean polymerization reaction of the diol and blowing reaction of water molecules. [OH] is the number of moles of the hydroxyl end groups per unit volume of polymerizing mixture and [W] is the number of moles of water molecules. The subscript 0 means initial value. The mass of each component per unit mass of the polymerizing mixture, ri, can be calculated from the concentration and the fractional conversion as follows.

MW and MCO2 are the molecular weight of water and carbon dioxide, respectively. rCO2,D is the initial mass of carbon dioxide that is dissolved in the polymerizing mixture, which is listed in Table 1. rBL can be obtained from the evaporation model for physical blowing agent, which will be discussed later. l is the heat of vaporization for blowing agent per unit mass, which is listed in Table 2 for CFC-11 (CCl3F) and HCFC-141b (CH3CCl2F). The energy balance equation states that the sum of the energy needed for temperature increase of the whole system and the heat consumption by the evaporation of the physical blowing agent is equal to the sum of the heat generated by the exothermic polyurethane reaction and the blowing reaction of the chemical blowing agent.

2.2. Gelling reaction

Reaction of diisocyanate and polyol ends up with gelling reaction of polyurethane, which is represented as below

If molecular diffusion is neglected and second order kinetics are assumed [8], the kinetic equation of gelling reaction can be represented as


AOH is the pre-exponential factor, EOH activation energy of the gelling reaction, Rg gas constant. The double underlined part of Eq. (8) denotes the effect of dilution of the reactive group due to the presence of the physical blowing agent and water in the reacting mixture [9]. The reaction parameters are listed in Table 3.

2.3. Blowing reaction

As water is used as the chemical blowing agent for polyurethane foam system, it reacts with isocyanate to form carbon dioxide and urea, which is the blowing reaction as shown below.


It has been generally accepted that the rate of water– isocyanate reaction is independent of the isocyanate concentration and follows the first order kinetics. However, Baser and Khakhar [9] reported that the second order kinetics predicted the water–isocyanate reaction more accurately than the first order kinetics. For this study, it is assumed that the reaction follows the second order kinetics as follows.


where the reaction parameters are listed in Table 3.

2.4. Evaporation of physical blowing agents

It is assumed that evaporation rate of the physical blowing agent is controlled by the heat generated due to chemical reaction. It is based on the fact that the rate of mass transfer is so fast that the blowing agent in the gas phase can be in equilibrium with the blowing agent in the liquid phase. Boiling point (TB) of the physical blowing agent in the reacting mixture depends on its mole fraction (xBL).
the mole fraction is a linear function of boiling temperature [8].
The relationship between mass ratio and mole fraction of the blowing agent in a liquid phase can be given as
where MB is the molecular weight of the blowing agent and Mno is the initial number average molecular weight of the polymerizing mixture listed in Table 1.
The material parameters of the CFC-11 and HCFC-141b for this model are listed in Table 2.

2.5. Density model

With the assumption of ideal gas, the density of the free rising foam at any time is given by

where p is the ambient pressure for the free rising foam. If the foam is pressurized during mold filling, the pressure should be changed to an appropriate value.

2.6. Viscosity model

An accurate viscosity model should be used for the analysis of polymer processing because the viscosity is one of the most important processing parameters. Castro–Macosko type model [15] is used in this study for numerical simulation of expanding foam.


Here, XNCO is isocyanate conversion and XNCO,g is its gel conversion. Parameters of the viscosity model are listed in Table 4. If the bubble size is small enough for the foam to be considered as a dilute bubble suspension, the foam is considered as a generalized Newtonian fluid as follows.

Ca, capillary number, is defined as (20)


where m is the viscosity of the matrix, G is the shear rate, R is the radius of the spherical bubble, and G is the interfacial tension. If the bubble size is uniform and the number of bubbles per unit polymer mass (nb) is known, the bubble radius can be obtained as
where Vg is the volume of gas phase per unit polymer mass and obtained as
In this study, nb is set to be 109 kgK1 according to the calculation by Niyogi et al. [18].

3. Numerical formulation

Prediction of mold filling by the self-expanding polyurethane foam was carried out based on the above theoretical modeling, such as energy balance, polyurethane reaction, blowing reaction, and evaporation of physical blowing agent. With the assumptions of ideal mixing and rapid bubble nucleation, the foam was modeled as a continuum. It was assumed that the continuum is a generalized Newtonian fluid whose constitutive equation is given by the foam rheology. For numerical calculations, a pressure based finite volume scheme was developed for unstructured meshes and the SIMPLE algorithm was employed with treatment of fluid compressibility. Cell based, co-located storage is used for all physical variables. In order to deal with the moving interface, an explicit high-resolution scheme that was similar to the CICSAM (compressive interface capturing scheme for arbitrary meshes) method was used. More detailed numerical schemes are described in the previous paper [14]. In addition, the multi-grid algorithm [17] was incorporated into the numerical code to increase the rate of convergence and reduce the calculation time compared with equivalent single grid schemes. The numerical code was written for fully three dimensional analysis by using the programming language, C++. The computational domain in the mold cavity is partially filled with initial charge of the foam and the moving interface is considered. It was assumed for the numerical simulation that the empty regions in the mold were filled completely with a fictitious fluid that had different physical properties from the foam. If surface tension at the interface between the foam and the air is negligible, the general governing equations for compressible Newtonian fluid with Stokes’ hypothesis include the continuity equation, momentum equation, and the energy equation as shown below

Here r is density, t is time, v is the velocity vector, g is the acceleration due to gravity, p is pressure, m is the shear viscosity, C is specific heat, and k is thermal conductivity.
Hg is heat generation during the polyurethane foam formation, which is the term on the right hand side of the Eq. (1). In addition to these equations, the mass conservation Eqs. (8), (9), and (15) are solved simultaneously. The time derivative v/vt in the mass conservation equations must be changed to the material derivative D/Dt. In order to track the interface between foam and air, the fractional volume function f(x,t) is defined such that

This function is governed by the following scalar advection equation.



The interface is located within the cells whose average value of f lies between 0 and 1. For these cells, material properties such as viscosity, density, specific heat, and thermal conductivity are interpolated linearly using the value of f.

For numerical calculations of the mold filling it is assumed that a part of the mold is initially filled with polymerizing mixture at the beginning of mold filling. The initial velocity is set to be zero at every point of the initial charge and the initial pressure distribution is uniform. The no-slip condition is applied at the mold boundary.

The overall numerical procedure for predicting mold filling of the polyurethane foam is as follows:

1. Fractional volume function f is initialized at all cells in the mold and the boundary conditions are defined.
2. All variables, such as pressure, velocity vector, temperature, density, viscosity, thermal conductivity, specific heat, and concentrations are initialized at all interior and boundary cells according to the fractional volume function.
3. Momentum equations and pressure correction equation are solved to obtain the velocity vector and pressure fields.
4. Mass conservation equations are solved to obtain the conversion of each species.
5. Energy equation is solved to obtain the temperature field.
6. Fractional volume equation is solved.
7. Density and viscosity of the foam are calculated by using the proposed models.
8. Density, viscosity, thermal conductivity, and specific heat are updated in the mold according to the newly calculated fractional volume function.
9. Time step is advanced and the numerical process is returned to step 3.
4. Results and discussion

4.1. Adiabatic stationary expansion

Self-expansion of polyurethane foam in an adiabatic cavity is studied numerically by neglecting momentum flow of the foam. Assuming that there is no resin flow or heat transfer through the mold boundary, governing equations given by the theoretical models are solved simultaneously by the numerical method. With the initial temperature and concentration of reactants and blowing agents, Eqs. (8), (9), and (15) are calculated and then Eq. (1) is treated by the explicit Euler scheme. Using the calculated


conversions and temperature variations, density and viscosity fields are predicted from Eqs. (16) and (17) at each time step.

Temperature variation of the polyurethane foam with respect to different blowing agents is shown in Fig. 1(a) when initial temperature is 50 °C. When 1.9 wt% of water is added to polyol as a chemical blowing agent, temperature increases more rapidly because the blowing reaction between water and isocyanate is exothermic. There is a synergy effect between the blowing reaction and the gelling reaction by sharing the generated heat. But when 24 wt% of CFC-11 is added to polyol as a physical blowing agent, temperature increase is retarded because the heat generated by the gelling reaction is consumed to evaporate the physical blowing agent. Therefore, the temperature increases moderately when both physical and chemical blowing agents are present. Fig. 1(b) shows the viscosity variation of the polyurethane foam with respect to the type of blowing agent when initial temperature is 50°C. The viscosity increases rapidly as the conversion of isocyanate reaches the gel point. The chemical blowing agent accelerates the polyurethane reaction by raising the temperature with the energy generated by the exothermic reaction. But the physical blowing agent decelerates the polyurethane reaction by absorbing the energy. Fig. 2 shows the density change of free rising foam with respect to different type of the blowing agent. Experimental data were obtained from the case that HCFC-141b and water were added by 24 and 1.9 wt%, respectively. Although the physical blowing agent used for the prediction is CFC-11, the density variation predicted from the theoretical model is compared with the experimental data. The experimentally measured density has the initial value of about 900 that is lower than the predicted initial density as shown in Fig. 2 because it takes some time for the mixed foam to be transferred from the mixing chamber to a cylinder in which the density is measured. The measured density becomes higher than the predicted value as time elapses because the gas generated from the blowing agent may diffuse to the outside of the foam in the experiment. Gas diffusion through the boundary of the foam is not considered for theoretical prediction.
4.2. Free rising in a thin block

The free rising of polyurethane foam has been studied by some researchers [8,11] in order to obtain information on the kinetics of the blowing and gelling reactions. The mold geometry employed in this study for numerical calculation of free rising polyurethane foam is hexahedral and the dimension is 0.02!0.2!0.1 m3, which is divided into 10,240 brick elements. Initially unexpanded foam is placed in the part of the mold, where z!0.02. The top surface at zZ0.1 is open and the other boundaries are solid walls


Initial temperature of the unexpanded foam is 50 ° C. The foam is blown by water which has been added to polyol by 1.9 wt% as a chemical blowing agent. It is assumed that the wall boundaries are adiabatic and the gravity force is acting in negative z direction. Fig. 3 shows the pressure contours and velocity vectors within the plane at yZ0.1 for free rising of polyurethane foam in the adiabatic cavity. The gravity effects on the shape of flow front are so significant that the gravity force can make the flow front flat. Before gelation of polyurethane, viscosity of the foam is about 1 Pa s. The maximum pressure at the corner is predicted as 1.03 Pa when the gravity is neglected, while the maximum pressure in the mold is calculated as 184.8 Pa when the gravity is taken into account. The pressure contours obtained by neglecting the gravity are similar to those of self-expanded fluid whose density decreases exponentially. Mold filling by a self-expanding fluid was studied previously by applying the pressure based finite volume method with the assumption that the foam density decreased exponentially with respect to time. Traces of fluid particles in the plane at yZ0.1 are shown in Fig. 4 for self-expansion of polyurethane foam in the adiabatic cavity. When the gravity force is neglected, the particles have the same traces as those predicted in the case of self-expanding foam whose density decreases exponentially. Every particle in the internal flow region moves to the central region, catches up with the flow front, and experiences the fountain flow. The fluid particle experiencing the fountain flow moves to regions near the wall and then slowly moves to the center region again. On the other hand, when the gravity force is acting in negative z direction, the particles change their paths. The particles near the bottom wall do not go up due to the gravity but just go to the center region. The particles in the other regions reach the flow front and then move to the wall more quickly than those in the case of without gravity.

4.3. Three dimensional filling of a cavity

Fig. 5 shows the mesh employed for numerical analysis of the mold filling by the expanding polyurethane foam. Unexpanded foam is placed initially within the space




defined by x!0.150. The boundary at xZ0.8 is an exit and the other boundaries are solid walls. Initial temperature of the unexpanded foam is 50 8C. The foam is blown by water, which has been added to polyol as a chemical blowing agent. The isothermal temperature boundary condition is assumed so that the wall temperature is maintained at 25 °C. For numerical analysis of the mold filling, 6273 nodes and 5120 finite volume elements are used. Flow front advancements are compared between the numerical and the experimental results for three different types of cavities as shown in Fig. 6. The figures on the left hand side show the numerical results in the plane at zZ0.02 for the water blown polyurethane foam with the isothermal wall when the gravity is present while those on the right hand side are the experimental results of the water and HCFC-141b blown polyurethane foam. The flow front shapes calculated from the numerical simulation are almost the same as those observed by the experiment, but the positions of the flow front are different. In the experiments, the polyurethane foam was mixed in the homogenizer and then the mixture was transferred to the mold cavity. Therefore, the foam was expanded slightly before it was transferred to the mold and the viscosity of the foam in the experiments was higher than that used in the numerical simulations. The observed flow fronts will move more slowly than the flow fronts predicted theoretically.

4.4. Filling of a small refrigerator cavity

As a practical problem, filling of a small refrigerator cavity between the ABS liner and the steel exterior panel is calculated numerically for the water blown polyurethane. Three dimensional shape of the refrigerator is depicted in Fig. 7. There are 48,185 nodal points and 38,144 cubic finite volume elements. Unexpanded foam is poured into the space defined by 0!x!0.24, 0!y!0.2, and 0!z!0.04. After it is supplied, the gate is closed. The exit is located in the plane of zZ0.4 at the upper part of the mold. It is assumed that the cavity has isothermal boundary and the gravity is acting in negative z direction. Fig. 8 shows the flow front advancements on the surface of the wall in two different directions. In the beginning of mold filling unexpanded foam flows out in x and y directions due to the gravity. After a while, the foam expands and fills up the cavity. It flows slowly at the corner and along the edge of the cavity so that the filling of the corner takes longer than that of other places. When the elapsed time is around 35 s, gelation of the foam is initiated and the viscosity increases rapidly. Fig. 9 shows the pressure contours and velocity vectors on the planes of xZ0.02, yZ0.02, and zZ0.02 when the elapsed time is 35 s. The highest pressure obtained is 450 Pa and located in the initially filled region. The speed of flow increases as the flow front approaches to the exit and the










distance between the pressure contours decreases near the flow front. The temperature contours on the planes of xZ0.02, yZ 0.02, and zZ0.02 are shown in Fig. 10(a) when the filling time is 35 s. The highest temperature is 510 K and observed in the initially filled region. The polyurethane foam is blown by water and the temperature increases significantly as mentioned previously. Density contours on the planes of xZ 0.02, yZ0.02, and zZ0.02 at 35 s are more complicated than the temperature and pressure contours as shown in Fig. 10(b). Foam density is dependent on temperature, pressure, and conversion of the blowing agent. In general, the foam density is inversely proportional to conversion and temperature, while it is proportional to pressure. The lowest density is about 50 kg/m3 and observed in the initially filled region. On the other hand the highest density is about 300 kg/m3 and occurs at the flow front. The highest density is about 6 times larger than the lowest value. According to Marciano et al., thermal conductivity of the CFC-11 blown polyurethane foam depends on the foam density as follows if the density is higher than 50 kg/m3.


With the relationship, the thermal conductivity varies from 0.016 to 0.045 W/m K, where the foam density is between 50 and 300 kg/m3. The thermal conductivity at the flow front is about 3 times larger than that at the initially filled region. Since, uniform temperature distribution is important for energy efficient control of the refrigerator, thermal resistance of the polyurethane foam in the thickness direction must be uniform throughout the entire cavity. Prediction of the cavity filling and foam density will be utilized for design of the refrigerator geometry and determination of the processing conditions.

5. Conclusions

A model was proposed for rigid polyurethane foam processing by assuming that evaporation of the physical blowing agent is controlled by heat generation and the carbon dioxide generation is controlled by the water–isocyanate reaction. To consider viscosity variation during foam expansion, a Castro–Macosko type model was adopted and foam rheology was applied to the model. Water, as a chemical blowing agent, makes the temperature increase more rapidly because blowing reaction between water and isocyanate is exothermic. However, the physical blowing agents retard the temperature increase because the heat generated from the gelling reaction is absorbed by evaporation of the physical blowing agent. The gravity effect on the flow front advancement and the pressure profile is significant because viscosity of the foam is very low before gelation. Three-dimensional filling of the cavity by the expanding foam was simulated numerically. The flow front advancement predicted by the numerical calculation agreed well with the experimental results. Finally, mold filling of water blown polyurethane foam was investigated numerically for a small refrigerator cavity. The density and thermal conductivity of the foam in the flow front were higher than those in the initially filled region.