# Numerical simulation of bubble growth in viscoelastic fluid with diffusion of dissolved foaming agent.

INTRODUCTIONIn recent years the quality control requirements for plastic foamed products have become increasingly stringent. For instance, for practical uses of microcellular plastics, supercritical fluid is required with a cell size on the order of microns, and a sufficient volume expansion ratio. Heat-resistant foamed trays made of polypropylene (PP) are required for many applications. Therefore, researchers have attempted to improve the foaming performance of polymer resins, and to quantify the rheological mechanisms of the foaming process.

Many aspects of bubble growth are strongly related to the performance of foamed polymer melts in the molding process. It is often necessary to control the bubble growth. For instance, it might be necessary to obtain foam in which small cells are uniformly distributed by suppressing the bubble growth, or to prevent a rapid volume expansion and coalescence of bubbles that cause the surface of the foam-molding goods to deteriorate. The bubble growth is induced by mass diffusion of dissolved gas, and is checked by the resistance of deformation in the liquid. However, many other intertwining factors are connected with this behavior, and thus it is difficult to specify the main factors that dominate the bubble growth.

Since the early days of plastic foamed products, many researchers have focused on the bubble growth with diffusion of dissolved gas in liquid [1, 2]. However, few papers have discussed the viscoelastic effect on the behavior of molten polymer. Venerus et al. [3] carried out a viscoelastic simulation of bubble growth with the Phan-Thien Tanner (PTT) model. They compared some approximation methods used for the simulation and discussed the accuracy of the results. They pointed out that the bubble growth rate is not significantly affected by the nonlinear viscoelastic characteristic. Agarwal [4] simulated bubble growth and collapse using the pom-pom constitutive equation. He then discussed the effect of the construction of branched polymers on the behavior. There are other noteworthy works related to bubble growth or collapse in viscoelastic fluid [5-9]. Those simulations have been done by using the approximation technique such as thin-boundary layer approximation around the bubble. Furthermore, they treated the liquid as infinite and did not express the saturation to the final bubble size during the later stages of vaporization. As a more practical simulation, diffusion-induced bubble growth models concerning the finite liquid surrounding the bubble have been developed [10-12]. Arefmanesh and Advani [13] and Ramesh et al. [14] carried out bubble growth simulations of viscoelastic fluid with a numerical model using the upperconvected Maxwell (UCM) model. Later, Lee and coworkers [15] performed a nonisothermal viscoelastic simulation of UCM fluid during the foam sheet formation, and discussed the importance of the cooling effect and the sheet thickness dependency. Ramesh and Malwitz [16] considered the concentration-dependent viscosity and diffusion coefficient for a nonisothermal bubble growth simulation [16]. Recently, Everitt et al. [17] extended this to a bubble growth simulation in reacting polymer using the Oldroyd-B model. To express the realistic bubble growth behavior in molten polymers for commercial use, viscoelasticity should be characterized in a more practical way as concerning multirelaxation modes within a wide time scale by the use of a viscoelastic constitutive equation that can also express the nonlinear property. Moreover, a limited amount of blowing agent dissolved in the molten polymer must be considered in order to express bubble growth in the later stages.

In this study, a new numerical model was developed that reflects these considerations [18]. A spherically symmetric finite element method (FEM) with Lagrangian nodes was used in the numerical model. The numerical scheme is advantageous for estimating the various distributions around the bubble, such as stress, strain, and dissolved gas concentration. Moreover, troublesome convection terms in the equations can be eliminated in the simulation procedure. The multi-relaxation mode PTT model was used to simulate bubble growth in molten PP. In the numerical study the effects of the rheological characteristics on the bubble growth behavior in the foam extrusion process of PP were investigated, as well as other controlling factors such as the amount of blowing agent, diffusion coefficient, surface tension, temperature, and ambient pressure.

Plastic foamed products have a nonuniform cell size distribution because of the nonuniform conditions and the frequency of the bubble nucleation in the polymer melt, which must be considered in predicting the cell size distribution. Significant studies about the numerical simulation of bubble growth have taken into account the distribution of nucleation in a nucleation model [19, 20]. In the current work we did with such a distribution related to nonuniform nucleation and the nucleus population density as the initial conditions for the bubble growth simulation. Rather, this investigation focused on specifying the main factors behind the bubble growth. Nevertheless, we studied how the nucleus population density influences the bubble growth and its volume expansion under the assumption of homogeneous cell nucleation.

NUMERICAL METHOD

Analysis Model

A numerical simulation model was constructed to illustrate a bubble in polymer melts after the nucleation generates consistently in the same interval and at the same time (Fig. 1). The radius of the spherically approximated region that influences the growth of a bubble is expressed as follows:

[S.sub.0] = (3/[4[pi]N])[.sup.1/3] (1)

[FIGURE 1 OMITTED]

where N is the nucleus population density. Equation 1 defines the size of a region before the expansion of the bubbles. A spherically symmetric, one-dimensional FEM simulation program was developed with following assumptions:

1. Because of the highly viscous entangled polymer melt, the inertial and gravity effects are ignored.

2. The effect of dissolved gas on the viscoelastic characteristic is ignored.

3. Before a bubble generates, there is a supersaturated dissolved gas state and the concentration distribution is uniform in the analyzed region.

4. Dissolved gas does not go in and out at the outer boundary of the analyzed region.

5. The polymer melt is incompressible. The volume of dissolved gas in the polymer melt is negligible.

6. The volume of gas in a bubble follows the ideal gas law.

7. The bubble contains only one component of gas.

8. Because the time scale of the cooling process in foaming is sufficiently larger than that of bubble growth, the gas in the bubble and the polymer melt surrounding it are regarded as being in an isothermal condition.

Governing Equations

Force Balance. Considering the viscoelastic stress, surface tension, pressure in the bubble, and the ambient pressure, the equation of force balance is expressed as

[[integral].sub.R.sup.S] [2/r] ([[tau].sub.[theta][theta]] - [[tau].sub.rr])dr + [P.sub.L] - [P.sub.G] + [[2[gamma]]/R] = 0 (2)

where r is the radial location; R is the radius of bubble; S is the radius of outer surface of the analyzed region; [[tau].sub.[theta][theta]] and [[tau].sub.rr] are the normal stresses in the tangential and radial directions, respectively; [P.sub.L] and [P.sub.G] represent the ambient pressure and the pressure inside the bubble, respectively; and [gamma] is the surface tension coefficient.

Diffusion of Dissolved Gas. Bubble growth is induced by the gas flux to the bubble, as a result of the diffusion of dissolved gas in the molten polymer. The mass transfer is expressed by

[[partial derivative]c]/[[partial derivative]t] = [[D.sub.G]/[r.sup.2]] [[partial derivative]/[[partial derivative]r]]([r.sup.2][[[partial derivative]c]/[[partial derivative]r]]) - [[[dot.R][R.sup.2]]/[r.sup.2]][[[partial derivative]c]/[[partial derivative]r]] (3)

where c, [D.sub.G], and t represent the dissolved gas concentration, diffusion coefficient, and time, respectively, and [dot.R] is the radial velocity at the interface of bubble to polymer.

Gas Concentration at the Bubble Interface. Generally, the relationship between pressure and gas concentration at the interface obeys Henry's law. This can be expressed as

[c.sub.R] = [k.sub.H][P.sub.G] (4)

where [c.sub.R] is the gas concentration at the bubble interface, and [k.sub.H] is Henry's constant.

Mass Balance of Gas. Under the assumption of the absence of gas flux at the outer interface, the total gas weight of the analyzed region in the outer surface is always constant. The mass balance of gas is expressed as

[[integral].sub.R.sup.S] 4[pi][r.sup.2]c dr + [[4[pi][R.sup.3][P.sub.G]]/[3[R.sub.GAS]T]] = [4/3] [pi][S.sub.0.sup.3][c.sub.0] (5)

where [c.sub.0] is the initial gas concentration in the analyzed region, [R.sub.GAS] is the gas constant, and T is the absolute temperature. The equation shows that the initial dissolved gas weight is equal to the sum of the gas weight in the bubble plus the gas weight in the polymer melt surrounding the bubble. This expression of mass balance is advantageous to avoid the treatment of gas concentration gradients at the bubble interface in the FEM simulation. In this study, the initial bubble radius [R.sub.0] was set at a very small value compared to the radius of the outer surface [S.sub.0]. The gas amount inside the initial bubble was negligible, and the total gas amount was calculated by the right side of Eq. 5.

Mass Balance of Polymer Melt. The volume change in the polymer melt by ambient pressure or gas dissolution is estimated to be sufficiently small in the present condition. It is recognized that these effects are not essentially important, considering the purpose of this study. To avoid complications, it was assumed that the polymer melt around a bubble is incompressible, and the volume of dissolved gas in the polymer melt is negligible. In other words, the melt density around a bubble is always constant. Therefore, a radial movement of an arbitrary position depends on the radial movement of the bubble interface, as expressed by

[r.sub.j.sup.3] - [r.sub.j0.sup.3] = [R.sup.3] - [R.sub.0.sup.3] (6)

where [r.sub.j] is an arbitrary location in the polymer melt, [r.sub.j0] is the initial location, and [R.sub.0] is the initial bubble radius.

Viscoelastic Constitutive Equation. To describe the rheological characteristics of the polymer surrounding a bubble, the multimode PTT model, which expresses the elongational behavior of polymer melts relatively well [21], was used in the numerical study. The equation is expressed as follows:

[Y.sub.i]([[tau].sub.i])[[tau].sub.i] + [[lambda].sub.i] [[[xi]/2] [[DELTA].[[tau].sub.i]] + (1 - [[xi]/2])[[nabla].[[tau].sub.i]]] = 2[[eta].sub.i]D (7)

[tau] = [n.summation over (i=1)] [[tau].sub.i]. (8)

For the function [Y.sub.i], two types were used, as expressed below:

PTT-I: [Y.sub.i]([[tau].sub.i]) = 1 + [[epsilon]/[G.sub.i]] tr ([[tau].sub.i]) (9)

PTT-II: [Y.sub.i]([[tau].sub.i]) = exp ([[epsilon]/[G.sub.i]]tr ([[tau].sub.i])) (10)

where [[tau].sub.i], [[lambda].sub.i], [[eta].sub.i], and [G.sub.i] indicate the stress tensor, relaxation time, zero shear viscosity, and relaxation modulus concerning relaxation mode i, respectively; [tau] is the total stress tensor; D is the rate of deformation tensor; [DELTA] and [nabla] represent lower- and upper-convective derivatives, respectively; and [epsilon] and [xi] are the model parameters of the PPT model. If [epsilon] = 0 and [xi] = 0, Eqs. 7-10 give the UCM model. In a spherically symmetric coordinate, Eqs. 6-10 for a particle whose initial location was at [r.sub.j0] can be expressed as follows:

[[[partial derivative][[tau].sub.[theta][theta]i]]/[[partial derivative]t]] = - [[Y.sub.i]/[[lambda].sub.i]] [[tau].sub.[theta][theta]i] + 2(1 - [xi]) [[[R.sup.2][dot.R]]/[[R.sup.3] + [r.sub.j0.sup.3]]] [[tau].sub.[theta][theta]i] + 2[G.sub.i] [[[R.sup.2][dot.R]]/[[R.sup.3] + [r.sub.j0.sup.3]]] (11)

[[[partial derivative][[tau].sub.rri]]/[[partial derivative]t]] = - [[Y.sub.i]/[[lambda].sub.i]][[tau].sub.rri] - 4(1 - [xi]) [[[R.sup.2][dot.R]]/[[R.sup.3] + [r.sub.j0.sup.3]]][[tau].sub.rri] - 4[G.sub.i] [[[R.sup.2][dot.R]]/[[R.sup.3] + [r.sub.j0.sup.3]]] (12)

[Y.sub.i] = 1 + [[[epsilon].sub.i]/[G.sub.i]] ([[tau].sub.rri] + 2[[tau].sub.[theta][theta]i]) or [Y.sub.i] = exp {[[[epsilon].sub.i]/[G.sub.i]]([[tau].sub.rri] + 2[[tau].sub.[theta][theta]i])} (13)

[[tau].sub.rr] = [n.summation over (i=1)] [[tau].sub.rri] (14)

[[tau].sub.[theta][theta]] = [n.summation over (i=1)] [[tau].sub.[theta][theta]i]. (15)

The strain-hardening property strongly depends on [epsilon], and the [xi] parameter is especially important if the second normal stress difference must be taken into account. In spherically symmetric bubble growth, the polymer melt around a bubble undergoes only biaxial stretching, and [epsilon] is enough to vary the nonlinear viscoelasticity. [xi] = 0 was set for all numerical runs.

Numerical Procedure

A one-dimensional FEM simulation for the radial direction was performed concerning the spherical symmetry of the analyzed object. Because the change in state would be sharper near the bubble interface, the division of FEM elements was constructed to become finer as the location got closer to the bubble interface. The number of elements was increased as the dependence of the numerical result on the division became negligible. The nodes on the polymer melt were set to move following their velocity. Displacements of these nodes were determined with Eq. 6. With the use of the Lagrangian method, the calculations of convection terms in the constitutive equation and the diffusion equation were eliminated. For the method of integration in Eqs. 2 and 5, the Gauss-Legendre integration was performed for each element and then summed up. The calculation of the diffusion of dissolved gas was done by the Galerkin method with linear interpolations in the elements. The natural boundary condition was set at the outer boundary of the analyzed region because of the absence of gas flux. Wilson's [theta] method was used for the time integration of the diffusion equation. The viscoelastic stress of each node was calculated by the integration of Eqs. 11-15 with the Runge-Kutta method in compliance with their rates of deformation. To integrate the constitutive equation, the time step size of stress calculation was often set to be smaller than that of the diffusion calculation, because the time step size must be sufficiently smaller than the relaxation time. The simultaneous equations described above were solved using the Picard iteration method for each time step.

Comparison With the Experimental Results

To estimate the validity of the numerical model, the numerical results of bubble growth of PP were compared with the experimental results obtained by Taki et al. [22]. The experimental procedure was performed as follows: A PP melt was attached to a high-pressure cell that was soaked in C[O.sub.2] gas under an isothermal condition (170[degrees]C) and under high pressure (11 MPa) until it reached solution equilibrium. Then the pressure in the cell was released. At the same time, the observation of the bubble generation and growth was started using a high-speed charge coupled devices (CCD) camera with a microscope. The experimental apparatus and procedure are fully described in the literature [22]. The measured material properties of the PP necessary for the simulation, such as the diffusion coefficient and Henry's constant, were described in detail by Areerat [23]. Figure 2 shows an example of the comparison between the experimental and numerical results. The viscoelasticity of the PP was characterized by the PTT-I model with six relaxation modes ranging from 0.001 s to 100 s [18]. The bubble growth phenomena differed in accordance with the bubble generation time. The experimental results in Fig. 2 are for the bubbles that can be seen at 5.5 s after the decompression starts. During that period, the pressure in the cell becomes about 7.6 MPa, while the rate of the pressure release is an almost constant 0.354 MPa/s value. It is indicated that the observed cross section of a bubble is proportional to the time during the early stage of bubble growth. Eventually, the pitch of the bubble growth becomes more gradual. At the early stage, there is enough gas around the bubble to encourage growth. However, at a later stage the gas supply to a bubble decreases because the growth of the neighboring bubbles and the diffusion of the gas to them are progressing similarly. The numerical results are shown in Fig. 2 for the various initial radii of the analyzed region. These radii correspond to their nucleation density. At the early stages the bubble growth is not influenced by the radius, and a good agreement is obtained between the numerical and experimental results. On the other hand, the initial radius considerably influences the bubble growth behavior at the later stages. When the initial radius is 0.75 mm, the numerical result agrees closely with the experimental data. Certain time lags were observed in the bubble nucleation [22]. The dissolved gas concentration will never regularly be distributed, as assumed in this numerical model. However, the bubble growth behavior could be approximately predicted not only the early stage, but also at the later stage in regard to the finite region and the limited amount of dissolved foaming agent around the bubble.

[FIGURE 2 OMITTED]

Numerical Study for Extrusion Foaming

Conditions. To investigate the dominating factors for bubble growth in extrusion foaming, numerical simulations with the present model were carried out. The basic conditions for the simulation are shown in Table 1. These material constants were set in reference to measured values at 190[degrees]C [23]. In this condition, the volume expansion ratio approaches 3.3 during the later equilibrium state. Considering the rapid pressure drop near the die exit, and the bubble growth under constant atmospheric pressure in the extrusion foaming process, the transition of the outer pressure for the simulations was set as shown in Fig. 3.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

Model Materials. Some linear PPs were used in this work. The storage modulus G' as a function of angular frequency for these materials is represented in Fig. 4. The experimental data were measured by RMS800 (Rheometric Scientific). Sets of the relaxation moduli characterized for these polymers are shown in Table 2. The calculated results of G' from these relaxation moduli are represented with lines in Fig. 4.

Effect of the Strain-Hardening Property. To estimate the effect of the strain-hardening property on the bubble growth rate, simulations were done for PP-A by varying the elongational property. Figure 5 indicates the effect of [epsilon] on the calculated transient biaxial elongational viscosity by the PTT-I model. The result of the PTT-II model with [epsilon] = 0.5 is represented. The PTT-II model can express the drop of the steady-state elongational viscosity with the increase of the strain rate. As a result, the PTT-II model causes considerably lower steady elongational viscosity than the PTT-I model at a high elongation rate. Figure 6 shows that the numerical result of the bubble radius growth corresponds to these very different aspects in the elongational property in Fig. 5. However, the difference in the calculated bubble growth property among them is very small. While the steady biaxial elongational viscosity with [epsilon] = 0.01 is about 100 times larger than that with [epsilon] = 0.5 in Fig. 5, there is little difference in the bubble growth rate in Fig. 6. Venerus et al. [3] carried out a nondimensional bubble growth simulation with the PTT model, and pointed out that the nonlinear viscoelastic property does not seriously affect the bubble growth. The above results suggest that this could be true for polymers used in conventional extrusion foaming processes.

[FIGURE 5 OMITTED]

Figure 7 shows changes in the dissolved foaming agent concentration, the rate of deformation, and the relative deformation gradient for the tangential direction in the bubble growth process. [D.sub.[theta][theta]] and [F.sub.[theta][theta].sup.-1] are tangential normal components of the rate of the deformation tensor and the relative deformation gradient tensor, respectively. These distributions were estimated for the radial direction from the numerical result of PP-A with [epsilon] = 0.5. The relative deformation gradient corresponds to the elongation ratio for the tangential direction. Just after nucleation, dissolved gas concentration decreases rapidly. At the same time the slope of the dissolved gas concentration near the bubble interface becomes remarkably sharper. As a result, a violent diffusion of gas into the bubble is generated. The inclination becomes more gradual as the growth of the bubble progresses. Regarding the rate of deformation, the value near the interface is largest just after the nucleation. Then the value at the interface decreases and the overall pitch for radial direction becomes gradual with the bubble growth. The relative deformation gradient is very different from the above aspect. A singularly large value appears at the interface. The peak value becomes larger as the bubble growth progresses. However, the value at a short distance from the interface stays lower even at the later stages of the bubble growth. This phenomenon is due to the large compression in the radial direction caused by the large tangential elongation near the interface of the bubble. The radial compression ratio increases in proportion to the square of the tangential elongation ratio. Therefore, the part that experiences the large deformation remains very narrow in the radial direction, even at the final stage. Figure 8 shows the normal stress difference around the bubble in the case of [epsilon] = 0.01 and [epsilon] = 0.5 with the PTT-I model. It is found that the effect of nonlinear viscoelastity appears only at the part that is extremely close to the interface. The value at most parts in the radial direction is reflected by the linear viscoelastic property as a consequence of the relatively small strain. The resistance for the bubble growth depends on the integration of the normal stress difference in the radial direction, as represented in Eq. 2. It is supposed that only the deference in the stress peak at the interface is insufficient to change the bubble growth rate. Accordingly, the strain-hardening property would not have an important effect on the bubble growth.

Effect of Linear Viscoelastic Property. The results described above imply that the linear viscoelastic property may be important for the bubble growth. Therefore, numerical simulations were carried out for some PPs whose linear viscoelasic characteristics differed from each other in response to their molecular weight and its distribution, as shown in Fig. 4 and Table 2. The results are shown in Fig. 9. The numerical conditions, except for the relaxation spectra, were the same as described in Table 1. The effect of the nonlinear property was not taken into account because of the previous results. The PTT-I model was used with [epsilon] = 0.5 in all cases. As a reference, the result from the diffusion-controlled condition is indicated in Fig. 9, in which no viscoelastic stress and no surface tension appear. It was found that the results of bubble growth behavior accompanied the linear viscoelasticity. It is remarkable that the bubble growth gets slower as G' at a higher frequency gets larger. Figure 10 shows the stress distributions of the PP-A at 0.02 s and 0.1 s after nucleation, involving the result for each relaxation mode. The figure indicates that the relaxation mode around 0.01 s influences the calculated total viscoelastic stress, which deters the bubble growth. The effect of the long relaxation time on bubble growth is found in the relationship between PP-D and PP-E in Fig. 9. During the later stages, the rate of deformation decreases and the dominant relaxation time over the resistance stress increases. Then the two lines of the two materials cross. However, the longer relaxation time does not greatly affect the whole bubble growth behavior, because the bubble growth has already reached its final stage by that time. The relaxation modulus at 0.01 s correlates with G' at 100 rad/s. Therefore, we tried to relate the numerical results for various PPs involving PP-A to PP-E with G' (100), which is the experimental result of G' at 100 rad/s. Figure 11 shows the result. The vertical axis [t.sub.1.5] indicates the time when the volume expansion ratio [S.sup.3]/[S.sub.0.sup.3] becomes 1.5 after the nucleation. The larger [t.sub.1.5] corresponds to the slower bubble growth. The obvious correlation between G' (100) and the bubble growth speed is shown in Fig. 11 for PP with C[O.sub.2] gas. The rate of the bubble growth will change according to the kind of polymer or foaming agent used. However, it is supposed that a similar result would be obtained in such a case by considering the dominating relaxation mode over the bubble growth rate. Figure 11 also represents the results for the cases in which the initial dissolved foaming agent concentrations were changed. The results indicate that the initial concentration has a remarkable effect on the rate of bubble growth. When the concentration is taken to be a higher value in order to get a foamed product with a higher volume expansion ratio, it is more difficult to restrain the bubble growth by viscoelastic stress.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

Analysis of the Main Factors. To determine the main factors for bubble growth among the many control factors involving the rheological properties, a multivariate analysis was done with 116 runs and different variables. The relaxation moduli of the above-mentioned PPs were used in these simulations with the PTT-I model. The calculated value of [t.sub.1.5] (i.e., the time when the volume expansion ratio reaches 1.5) after the nucleation was selected as the criterion variable. This would provide useful insights into the apparent volume expansion behavior of foam in the actual foaming process. The results of the standard regression coefficient are represented in Table 3. In the case of varying temperatures, Arrhenius time-temperature shift factors were used with an activation energy of 9 kcal/molK. The ranges of variables were established respecting the ranges of measured values and the actual foaming conditions. A larger positive value of the standard regression coefficient in Table 3 means that the explanatory variable has a stronger positive correlation with the delay of the volume expansion of the foam. In contrast, a large negative value means that the volume expansion rate increases with the increase of the explanatory variable value. Some influential factors were found for the rate of volume expansion other than the storage modulus at 100 rad/s. These factors are the initial outer radius, the initial dissolved agent concentration, the ambient pressure, and the diffusion coefficient, followed by Henry's constant and temperature. We found that [epsilon] and surface tension had little effect on the behavior. Figure 12 shows the effect of the initial outer radius [S.sub.0] on the bubble growth and the volume expansion. The numerical conditions, except for [S.sub.0], were set as in Table 1. In the case of a smaller [S.sub.0], R reaches its equilibrium value earlier. That is due to the smaller amount of dissolved foaming agent in the analyzed region and the shorter diffusion distance around the bubble. On the other hand, the smaller [S.sub.0] corresponds to the larger number of nucleus. As the result, the volume expansion ratio grows more rapidly with the smaller [S.sub.0]. The case with the initial outer radius of 0.05mm corresponds to the ultimate cell radius of 0.0637 mm and the ultimate cell density of 548/[mm.sup.3]. On the other side, the case with the initial outer radius of 0.2 mm corresponds to the ultimate cell radius of 0.263 mm and the ultimate cell density of 9.12/[mm.sup.3]. The effect of the ambient pressure is represented in Fig. 13. In this case, the initial dissolved foaming agent concentration was set to 160 mol/[m.sup.3]. The other conditions, except for the pressure, were set as in Table 1. Figure 13 shows the serious effect of the ambient pressure on the bubble growth behavior. The higher ambient pressure slows the bubble growth and makes the radius smaller. In the case of the higher ambient pressure, the gas amount in the bubble decreases with decreasing a supersaturated dissolved foaming agent. Also the bubble becomes more compressed because of its high compressibility. Figure 14 shows the effect of the diffusion coefficient. The numerical conditions, except for the diffusion coefficient, are shown in Table 1. The diffusion coefficient also has an important effect on the bubble growth. It is naturally assumed by the fact that the origin of the bubble growth is diffusion. We have already pointed out the serious effect of the initial dissolved foaming agent concentration [c.sub.0] on the behavior in Fig. 11. The higher concentrations lead to more violent diffusion around the bubble. Figure 15 shows the effect of Henry's constant on the bubble growth behavior. The numerical conditions, except for Henry's constant, are the values in Table 1. A higher Henry's constant causes a somewhat slower bubble growth, as the result of the decreasing supersaturated foaming agent. In an actual foaming process, it is supposed that Henry's constant has an important influence on the controlling [c.sub.0] value in the mixing and dissolving process before the bubble is generated. In the analysis, the initial dissolved foaming agent concentration [c.sub.0] was set as the initial condition. In the case of a constant [c.sub.0], Henry's constant does not have much effect on the bubble growth, as shown in Fig. 15. Some important factors in the bubble growth and volume expansion were specified above. To interpret the bubble growth behavior in an actual foaming process, these main factors must be considered in addition to the rheological properties.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

CONCLUSIONS

Viscoelastic simulations of bubble growth in PP physical foaming under isothermal and constant ambient pressure conditions were performed. A multimode PTT model was used to analyze the dynamic growth behavior of spherically symmetric bubbles with the diffusion of a foaming agent (C[O.sub.2]). A change in the dissolved foaming agent concentration in polymer and a change in the strain of the polymer melt surrounding the bubbles were simulated with the Lagrangian FEM method. The simulation technique was validated by comparing the numerical results with the bubble growth rates, which were experimentally obtained by visual observation of the PP/C[O.sub.2] batch foaming system. Actual bubble growth behavior was expressed during both the the early and later stages, respecting the finite region and the approximate limited amount of dissolved foaming agent around the bubble. The simulation results clearly showed that the strain-hardening characteristic of polymer does not strongly affect the bubble growth rate. Though the strain-hardening property has important affect on the nonuniform deformation around the bubble when the bubbles approach mutually, it is not so important for the uniform deformation around the bubble. The linear viscoelastic characteristic is more influential, and the relaxation mode around 0.01 s is the important factor in determining the bubble growth rate in the early stage of foaming. A multivariate analysis was done to determine the main factors that contribute to the bubble growth stage among the many control factors involving the rheological properties. It was clarified that the bubble nucleus population density, the surrounding pressure, initial dissolved foaming agent concentration, and diffusion coefficient are more important factors than the viscoelastic characteristics.

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

ACKNOWLEDGMENTS

The authors thank Prof. M. Ohshima, Dr. S. Kihara, and K. Taki of the Department of Chemical Engineering at Kyoto University for providing the visual observation data regarding bubble growth, and for useful discussions on the subject.

REFERENCES

1. E.J. Barlow and W.E. Langlois, IBM J. Res Dev., 6, 329 (1962).

2. J.R. Street, A.L. Fricke, and L.P, Reiss, Ind. Eng. Chem. Fluid., 10, 54 (1971).

3. D.C. Venerus, N. Yala, and B. Bernstein, J. Non-Newtonian Fluid Mech., 75, 55 (1998).

4. U.S. Agarwal, e-Polymers, 14, 15 (2002).

5. E. Zena and L.G. Leal, Ind. Eng. Chem. Fund., 14, 175 (1975).

6. C.D. Han and H.J. Yoo, Polym. Eng. Sci., 21, 518 (1981).

7. C.D. Han and H.J. Yoo, AIChE J., 28, 1002 (1982).

8. R.K. Upadhyay, Adv. Polym. Technol., 5, 55 (1985)

9. J.S. Vrentas and C.M. Vrentas, J. Appl. Polym. Sci., 67, 2093 (1998).

10. M. Amon and C.D. Denson, Polym. Eng. Sci., 24, 1026 (1984).

11. A. Arefmanesh, S.G. Advani, and E.E. Mishaelides, Int. J. Heat Mass Transfer, 35, 1711 (1992).

12. D.C. Venerus, Polym. Eng. Sci., 41, 1390 (2001).

13. A. Arefmanesh and S.G. Advani, Rheol. Acta, 30, 274 (1991).

14. N.S. Ramesh, D.H. Rasmussen, and G.A. Campbell, Polym. Eng. Sci., 31, 1657 (1991).

15. S.T. Lee, N.S. Ramesh, and G.A. Campbell, Polym. Eng. Sci., 36, 2477 (1996).

16. N.S. Ramesh and N. Malwitz, J. Cell. Plast., 35, 199 (1999).

17. S.L. Everitt, O.G. Harlen, H.J. Wilson, and D.J. Read, J. Non-Newt. Fluid Mech., 114, 83 (2003).

18. Y. Otsuki, T, Kanai, K. Taki, and M. Ohshima, Seikei-Kakou (J. JSPP), 15, 638 (2003).

19. M.A. Shali, K. Joshi, and R.W. Flumerfelt, Chem. Eng. Sci., 52, 635 (1997).

20. Y. Inamori, K. Hayama, M. Takada, M. Ohshima, and M. Tanigaki, Seikei-Kakou (J. JSPP), 11, 194 (1999).

21. N. Phan-Thien and R.I. Tanner, J. Non-Newt. Fluid Mech., 2, 35 (1977).

22. K. Taki, T. Nakayama, T. Yatsuzuka, and M. Ohshima, J. Cell. Plastics, 39, 155 (2003).

23. S. Areerat, "Solubility, Diffusion Coefficient and Viscosity in Polymer/C[O.sub.2] Systems," Ph.D. Dissertation, Kyoto University, Kyoto, Japan (2002).

Yasuhiko Otsuki, Toshitaka Kanai

R & D Division, Prime Polymer Co., Ltd., 580-30, Nagaura, Sodegaura City, Chiba Pref., 299-0265, Japan

Correspondence to: Y. Otsuki; e-mail: yasuhiko.otsuki@primepolymer.co.jp

TABLE 1. Input data for bubble growth simulation in extrusion. Radius of analysis region 0.1 mm Initial radius of bubble 0.1 [micro]m Material constants Surface tension coefficient 1.32 X [10.sup.2] N/m Henry's constant 1.1 X [10.sup.-4] mol/Nm Diffusion coefficient 8.0 X [10.sup.-9] [m.sup.2]/s Condition Initial density of C[O.sub.2] 7.0 X [10.sup.1] mol/[m.sup.3] Temperature 190[degrees]C Atmosphere Pressure 1.0 X [10.sup.-1] MPa TABLE 2. Relaxation modulus used for numerical analysis (190[degrees]C). G (Pa) [lambda] (s) PP-A PP-B PP-C PP-D PP-E 0.001 60000 26000 35000 24000 120000 0.01 30000 21000 25000 17000 100000 0.1 13000 9500 15000 13000 45000 1 4100 3000 8800 13000 15000 10 750 480 4100 9000 2000 100 80 1 500 3000 230 1000 5 0 0 300 0 TABLE 3. Standard regression coefficients for the numerical result of [t.sub.1,5]. Standard regression Range coefficient Initial outer radius 5 X [10.sup.5]-2 X [10.sup.4] m 0.421 [S.sub.0] Initial dissolved 40-160 mol/[m.sup.3] -0.409 agent concentration [c.sub.0] Atmospharic pressure 0.1-0.5 MPa 0.401 Diffusion coefficient 2 X [10.sup.-9]-2 X [10.sup.-8] -0.286 [m.sup.2]/s Storage modulus at 0.021-0.128 MPa 0.227 100 rad/s Henry's constant 1 X [10.sup.5]-1.5 X [10.sup.4] 0.110 mol/Nm Temperature 160-220[degrees]C -0.060 Non-linear parameter 0.01-0.9 -0.006 [epsilon] Surface tension 5 X [10.sup.-3]-3 X [10.sup.-2] N/m 0.004

Printer friendly Cite/link Email Feedback | |

Author: | Otsuki, Yasuhiko; Kanai, Toshitaka |
---|---|

Publication: | Polymer Engineering and Science |

Date: | Sep 1, 2005 |

Words: | 5997 |

Previous Article: | Determination of the compatibility of NBR-EPDM blends by an ultrasonic technique, modulated DSC, dynamic mechanical analysis, and atomic force... |

Next Article: | Mechanical behavior of high impact polystyrene based on SBR copolymers: Part I. |

Topics: |