Effects of the Polymer Fiber on the Flow Field from a Slot Melt Blowing

Exxon slot melt blowing dies consist of dual, rectangular, converging jets and are used in the industrial melt blowing process to attenuate molten pol...
0 downloads 0 Views 232KB Size
Download PDF
Ind. Eng. Chem. Res. 2008, 47, 935-945

935

Effects of the Polymer Fiber on the Flow Field from a Slot Melt Blowing Die Holly M. Krutka, Robert L. Shambaugh,* and Dimitrios V. Papavassiliou School of Chemical, Biological and Materials Engineering, The UniVersity of Oklahoma, 100 East Boyd Street, SEC T335 Norman, Oklahoma 73019

Exxon slot melt blowing dies consist of dual, rectangular, converging jets and are used in the industrial melt blowing process to attenuate molten polymer fibers. The air flow creates a drag force that accelerates the polymer and rapidly reduces the fiber diameter. For previous experimental and computational fluid dynamics studies, the effect of the fiber on the air was assumed to be negligible. By including the fiber as a boundary in the computational domain, this assumption was tested. A modified version of the Reynolds stress model was used to simulate the turbulent air flow field. It was determined that the centerline air velocity (located halfway between the fibers) has an increased maximum due to the presence of the fiber. In addition, the turbulence in the flow field is dampened by the presence of the fiber. The jet spreading rate is higher halfway between the fibers than at the center of the domain, where the fiber is located. The air flow around the fiber is nonuniform, leading to varied shear stress at different radial positions on the fiber edge. The temperature of the air flow is related to the ratio of the polymer to air flow rate; as this ratio increases, the hot fiber reduces the cooling rate of the air. 1. Introduction Exxon slot dies are commonly used for the commercial production of fibers in the melt blowing process. This type of melt blowing die consists of two converging slot air jets that are located on either side of a row of polymer capillaries.1-3 With 20 or more capillaries per centimeter of die width, a typical 1 m Exxon die simultaneously produces thousands of filaments. These filaments are collected on a screen to form a nonwoven mat of fibers.3 Figure 1 shows an Exxon slot melt blowing die with a partial cutout.4 Figure 1 also shows a cross-section of the same slot melt blowing die. In the melt blowing process, molten polymer is forced through the capillaries by an extruder. Then, the air jets (shown in Figure 1) create a high drag force on the surface of the fiber that leads to rapid attenuation of the fiber diameter. The diameter of the fiber is often reduced by a factor of 1000 in a matter of microseconds.4 Due to their importance in numerous applications, planar jets have been researched in depth. Statistically planar, or rectangular, jets are two-dimensional. One well-documented characteristic of planar jets is their self-similar structure. This structure is observed beyond 30 slot widths away from the jet.5 For selfsimilarity, the velocity should be nondimensionalized using the nominal discharge jet velocity, Vjo. Also, the position must be nondimensionalized using x1/2, which is the distance from the centerline where the velocity is equal to half of the centerline velocity.6 The concept of self-similarity allows for the velocity profiles to be predicted. The centerline decay is related to z-1/2, where z is the primary direction of flow. Experimental verification of self-similarity is contained in the work of both Bradbury5 and Heskestad.7 They both used hot-wire anemometry to measure the velocity and turbulence profiles from a turbulent planar jet. For the air-side of a melt blowing die, previous work has shown that computational fluid dynamics (CFD) can be used to successfully match experimental measurements of both velocity and temperature profiles for slot dies8-11 and for annular dies.12-14 Research has also been focused on the understanding * To whom correspondence should be addressed. Tel.: (405) 3256070. Fax: (405) 325-5813. E-mail: [email protected].

of the fiber (polymer) side of melt blowing. For melt blowing, theoretical models have been developed to predict the fiber diameter, fiber temperature, and fiber stress as functions of position below the die. (The more advanced models can also predict fiber vibration.) This modeling has been done for both solid15-17 and hollow18 fibers. Uyttendaele and Shambaugh15 first developed a model in 1D; they found that using either a Newtonian viscosity or a Phan-Thien viscosity resulted in very similar predictions of final fiber diameters. Their results compared well with experiments. Rao and Shambaugh16 extended the model to include 2D; this extension allowed prediction of vibrations in the melt blowing process. Finally, in a further model advancement, Marla and Shambaugh17 extended the model to include 3D (rather than planar) fiber movement during the melt blowing process. The air velocity and temperature profiles from the simulations provided by Krutka et al.8-10 were input into the model as boundary conditions for the air-polymer interface. Other inputs to the model included air temperature, air flow rate, polymer temperature, and polymer flow rate. The model predicted the fiber temperature, velocity, diameter, frequency of vibration, and amplitude of vibration.17 Yin et al.19 combined on-line measurements with off-line analysis of web structure to provide a fundamental description of the melt blowing process. For a blunt slot melt blowing die, Wang and Ke20,21 used a model of spheres and weightless springs to simulate the movement of a fiber during melt blowing. For the simpler case of melt spinning (where the final fiber velocity is determined by a take up roll), Hietel and Marheineke22 used modeling and numerical simulation of fiber dynamics to predict air and fiber interactions. The interactions between the air and the fiber as they occur in melt blowing are difficult to model. In most previous work,15-18 the air velocity and temperature fields were entered as boundary conditions in the solution of the polymer side equations. The velocity fields were from experimental correlations that were independently measured23 or predicted by CFD.8,9 Similarly, the temperature fields were from experimental correlations24,25 or from CFD.10 For the experimental measurements of these temperature and velocity fields, no fiber was

10.1021/ie070871i CCC: $40.75 © 2008 American Chemical Society Published on Web 01/01/2008

936

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008

Figure 1. (a) “Exxon” slot melt blowing die (with a partial cutout) showing the row of polymer capillaries and the two air slots and (b) cross-sectional view of the slot die. For the die simulated herein the values of the parameters are Dc ) 0.475 mm, h ) 3.32 mm, b ) 0.65 mm, and θ ) 60°. Table 1. Polymer Flow Rates and Initial Air Velocities Used in the Simulations case

mf (g/min)

Vf,die (m/s)

Qair (slpm)

ma (g/min)

Va,die (m/s)

momentum flux ratio

Re

1 2 3 4 5 6 7

0 0.55 0.55 0.275 0.55 1.1 0.55

0 0.058 0.058 0.0289 0.058 0.1156 0.058

100 300 200 100 100 100 50

120 360 240 120 120 120 60

33.7 101.1 67.4 33.7 33.7 33.7 16.85

N/A 4.87 × 10-4 1.09 × 10-3 1.09 × 10-3 4.38 × 10-3 1.74 × 10-2 1.74 × 10-2

800 2300 1500 800 800 800 400

present during the measurements (i.e., the effect of the fiber on the field was not included). Likewise, for the CFD simulations, no fiber was placed as a boundary condition in the simulations. In contrast to the aforementioned work, Krutka et al.26 used CFD to study the effect of the presence of a fiber in an axisymmetric simulation of a single annular jet. In the present paper, fiber-inclusive simulations were conducted for another type of die geometry: the Exxon slot melt blowing die. The goals of the present paper are to (a) examine how the air flow field changes due to the presence of a polymer fiber, (b) determine how different air and polymer flow rates affect the flow and temperature field, and (c) examine the stress on the fiber due to the air drag.

condition for the CFD simulation. The velocity and temperature of the fiber, which vary with position, were specified for the CFD simulations with user defined functions (UDF). All the simulations discussed in this paper were run for the die type shown in Figure 1. The jet width at the die face, b, was 0.65 mm. The distance between the outside edges of both jets, h, was 3.32 mm. For all simulations, the angle θ between the die face and the jets was 60°. A melt blowing die with this configuration has been studied experimentally by Harpham and Shambaugh.23,24 Seven different simulations are discussed in this paper. The difference between the cases is in the air and/or polymer flow rates. The momentum flux ratio, which is used to distinguish the simulation cases, is defined as follows:29

2. Numerical Modeling and Simulation Parameters The computational fluid dynamics software package used to complete all the simulations discussed in this paper was Fluent 6.2.27 The determination of the geometry of the computational domain, the generation of the computational grid, and the turbulence modeling are discussed in detail in the following sections. 2.1. Problem Configuration. For simulations of an Exxon slot die without the presence of a fiber, the geometry is statistically 2D.8-11 However, the inclusion of a fiber requires a 3D simulation. The computational domain and grid were created in Gambit.28 To determine the diameter, velocity, and temperature of the fiber as a function of distance from the die, the model developed by Marla and Shambaugh17 was used. Then, the predicted fiber diameter was used as a boundary

qj )

Ff,dieVf,die2 Fa,dieVjo2

(1)

All cases are listed in Table 1. The case number increases as the momentum flux ratio increases. Case 1 is the only simulation that does not include a fiber (mf ) 0), and case 1 is used to compare the results of an air-only simulation to those simulations that include a fiber. Cases 1 and 4-6 have the same air flow rates as experiments completed by Harpham and Shambaugh23 (their experiments were completed without the presence of a fiber). 2.2. Grid Generation. Figure 2 shows the computational domain for all simulations. This domain extends 50 mm below the die face, in the z-direction. The domain also includes the

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 937

Figure 2. Computational domain used in the simulations and a close up of the area of interest in the computational domain. Boundary conditions are included in the figure. For the die simulated herein the domain width in the y-direction was 0.7257 mm.

air jet, which extends 5 mm above the die face. The width of the domain is 15 mm in the x-direction. The y ) 0 plane is at the center of the domain. The width of the domain in the y-direction is 0.7257 mm. This domain size allows the modeling of one polymer hole in a die with 35 polymer holes/in. (1380 holes/m). Although the simulation cannot be completed in 2D (due to the presence of the fiber), using symmetry for the boundary condition on the left side of the domain (x ) 0) allows the computational domain to be reduced. In previous CFD simulations, the grid was refined in Fluent in order to reduce the size of the cells in the area of most interest.8-11 This approach was also attempted with the present geometry, but a sharp change in the cell size close to the fiber created discontinuities in the simulation results, as well as instability during the simulations. It is desirable to have smaller cells in the area of most interest, which is the area close to the die face and close to the fiber, without rapid changes in the cell size. To accomplish this, the average cell spacing was specified on each edge in the domain with a successive ratio (the fraction of cell edge length increased) between consecutive cells. Next, the faces and domain volume were meshed using quadrilateral cells. On the right side of Figure 2 there is a close up of the section of the computational domain that includes the jet, fiber, and die face. Periodic boundary conditions were used on the front and back of the domain, in order to simulate the flow around polymer fibers in the center of the die (where die end effects are negligible). The fiber (with a diameter that changes with z-position) was represented using a round half-cylinder wall with a no-slip surface. The bilateral symmetry of the die about the x ) 0 plane was used to cut the domain in half (i.e., only one of the slot jets was considered). The boundary condition for the air entrance into the jet was a mass inlet. Because the simulation was non-isothermal and several different air flow rates were simulated, compressibility could not be neglected, and the mass inlet is the appropriate choice for this boundary (as opposed to

Figure 3. Centerline velocity at y ) Lc/2 for different grid resolutions. These simulations were all run for Vjo ) 33.7 m/s and mf ) 0.55 g/min. Table 2. Simulation Grid and the Computational Requirements for the Simulations Shown in Figure 2 no. of grid cells

no. of iterations

computation time (h)

653 090 1 023 162 2 640 102

15 177 18 736 48 000

15 29 195

a velocity inlet boundary that was used for 2D isothermal simulations8). The die face was represented as a no-slip wall. Several different simulations were run at different grid refinements in order to determine the cell size necessary to achieve grid independence. Figure 3 shows the dimensionless velocity at the y ) Lc/2, x ) 0, centerline for three different grids. Table 2 gives the number of cells, the number of iterations necessary, and the computational time when the simulation was run using four Dell Pentium 4 Xeon64 processors in parallel. (Note that the simulation using the finest grid refinement did

938

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008

not converge to 10-5 as did the other two simulations; instead, this finest grid simulation was stopped after it achieved 5 × 10-5 convergence.) For the grid with 653 090 cells, the symmetry (left boundary) and side pressure outlet (right boundary) were split up into five sections. These sections had increasingly larger cells with increasing z-position. The quadrilateral cells close to the die face were the smallest; these cells had an edge length of 0.10 mm. The largest cells, at the bottom of the domain, had an edge length of 0.22 mm. The cell width in the y direction was 0.2 mm. For the 1 023 162 cell simulation, the edges of the computational domain were specified using an average cell size of 0.08 mm with a successive ratio of 1.002. The cells closest to the die face had an edge length of 0.04 mm, while the cells at the bottom of the domain (far from the areas of interest and far from high gradients) had cells with a side length of 0.14 mm. The cell width in the y direction was 0.12 mm. Although this is a wide range of cell sizes, the change in the cell size is gradual, which is important for the simulation results and stability. Finally, the finest grid had 2 640 102 quadrilateral cells. This grid contained cells with an edge length of about 0.03 mm for locations close to the die face. Farther from the die face, the edge length was 0.14 mm and the cell width in the y direction was equal to 0.1 mm. Figure 3 shows that all three simulations resulted in a centerline velocity with a maximum close to z/h ∼ 1. However, the maximum reached for the case with 653 090 cells was slightly lower than the other two cases. Therefore, the grid with 1 023 162 was used for all simulations. Figure 2 shows this computational domain and the edges labeled according to the cell size. Sides with the label “A” have an average cell edge length of 0.08 mm with a successive ratio of 1.002. The walls and inlet of the air jets are labeled “B”; these sides have cells with 0.08 mm spacing. The sides labeled “C” have cells with 0.12 mm spacing. These cell sizes were used for all the simulations listed in Table 1. 2.3. Turbulence Modeling. Krutka et al.10 have shown that the ideal gas equation can accurately represent the compressibility in non-isothermal simulations and simulations with high air flow rates similar to those discussed in this paper. Therefore, for all simulations listed in Table 1, the ideal gas equation was used to model the air density. All simulations were run with an inlet air temperature of 330 °C (the inlet excess temperature was Θjo ) 309 °C). The inlet polymer temperature was 295 °C. These temperatures were based on the experiments by Harpham and Shambaugh23 and the modeling work of Marla and Shambaugh.17 Krutka et al. have shown that the Reynolds stress model (RSM) can accurately simulate the flow field from different configurations of slot melt blowing dies.8-11 Several other Reynolds averaged Navier-Stokes (RANS) equation models have been compared to experimental results (including standard k- and realizable k-), and only the RSM-obtained results have been validated with experiments for use in the case of plane slot dies.8-11 Since the presence of the fiber is expected to disturb the symmetry of the jet flow and produce anisotropic intensities of the turbulent velocity fluctuations, modeling the individual components of the Reynolds stress separately, as is done in the RSM model, is expected to be more accurate than using other models such as the k- model, where the components of the Reynolds stress are grouped together under the turbulent kinetic energy term. Therefore, the RSM was used, with the modifications suggested by Krutka et al.8-11 to simulate the air flow for all the simulations in Table 1. The following equation models the transport of Reynolds stresses:27

∂ ∂ (FUkuiuj) ) - [Fuiujuk + p(δkjui + δikuj)] + ∂xk ∂xk ∂Ui ∂Ui ∂ ∂ µ (u u ) - F uiuk + ujuk + ∂xk ∂xk i j ∂xk ∂xk

[

] (

p

(

)

)

∂ui ∂uj ∂ui ∂uj + - 2µ (2) ∂xj ∂xi ∂xk ∂xk

Equation 1 is written using the summation convention. The rate of dissipation is modeled using the following equation:27

F

D ∂ ) Dt ∂xi

[( ) ] µ+

µt ∂  2 + C1 Gk - C2F σ ∂xi k k

(3)

The Fluent software provides default values for C1 and C2 of 1.44 and 1.92, respectively.27 However, previous studies have determined that using these parameters does not lead to an accurate reproduction of experimentally measured flow fields. Using the default values of these constants results in lower mean air velocity close to the die face and higher mean velocity away from the die face. To reduce the dissipation of the turbulent kinetic energy and, thus, to increase the mean velocity, the second term of the right-hand side of eq 3 has to become smaller and the third term has to become larger (in absolute value). Therefore, the value of C1 was decreased to 1.24, and the value of C2 was increased to 2.05 (these modified values were determined by comparing experimental data with the simulations results8). The enhanced wall treatment option was enabled for the simulations in Table 1. To more accurately model the conditions during melt blowing, the simulations were completed for non-isothermal conditions. Enabling the energy equation was necessary to simulate the air temperature field. The energy equation used in the simulations was as follows:27

∂ ∂ (FE) + [U (FE + p)] ) ∂t ∂xi i ∂ ∂xi

[(

kc +

)

]

cpµt ∂T + Uj (τij)eff (4) Prt ∂xi

The density of the gas was calculated using the ideal gas law. To decrease the amount of computational time required and increase the stability of the simulations, isothermal conditions with constant air density were applied until the residuals of the model equations reached 10-4. Then, the density was determined using the ideal gas equation, and the appropriate temperature was applied at the boundaries. After these changes, the simulations were run until the residuals for velocity reached the value of 10-5 convergence, except for the energy equation, which was run until the temperature residual values reached 10-6. The residuals were defined as the root mean square of the residuals in each cell of the domain, normalized by dividing the residual value at every iteration by the maximum residual value that occurred during the first five iterations. 3. Results and Discussion 3.1. Results from the Fiber Model. For cases 2-7 the diameter of the polymer capillary was 475 µm. This dimension was used as the diameter of the polymer at the die face, and polymer diameter was assumed to decrease for z-positions below this point (die swell was assumed small). The radius as a function of z-position was calculated using the model developed by Marla and Shambaugh17 (see Figure 4). As expected, the

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 939

Figure 4. Fiber radius as predicted with the rheological (polymer side) model.

Figure 5. z-velocity as predicted with the rheological model.

higher the momentum flux ratio, the larger the fiber diameter. In addition, the model predicts, within the range of values tested, that the momentum flux ratio, not the air flow rate, determines the fiber diameter. This is based on the fact that cases 3 and 4, which have the same momentum flux ratio, but different air and polymer flow rates, have similar diameters. In addition, cases 6 and 7, which also have the same momentum flux ratio, have similar diameters. Figure 5 shows the velocity of the fibers, as predicted by the Marla and Shambaugh17 model. The fibers with smaller diameters are traveling at higher velocities. The size of the computational domain used in the CFD simulations was based on these radii. Specifically, for a particular air flow rate, the radius profile from the polymer side model was used to construct the domain for the CFD simulation at the same air flow rate. In addition, the fiber velocities were read into the simulation through user defined functions and the wall representing the fiber was set as a moving wall. Since the simulations were not time-dependent, the mean flow of the air and polymer was modeled. This approximation ignores any fiber oscillation caused by air turbulence. 3.2. Effects of the Polymer on the Mean Flow Field. Without a fiber present, dual rectangular jets have been found to exhibit three regions of development.30 In the first region the jets are separate and a maximum velocity is present due to each individual jet. High mean velocity gradients and Reynolds

Figure 6. Location of the maximum air velocity at different z-positions for cases 1 and 4-6 on the y ) 0 plane.

stresses between the two jets lead to jet interactions. For the blunt melt blowing die, a recirculation area is present between the two jets.8 This recirculation area results in a negative centerline air velocity. The second region in the flow field is characterized by the beginning of the jets merging together. The centerline velocity is positive in this region, but there are still two separate velocity maxima (one due to each jet); these maxima are closer to each other than in the first region. The second region exhibits a maximum in the turbulence intensity and turbulence kinetic energy. The flow in the third region can be described as self-similar. In this region, there is only one velocity maximum, and this maximum occurs on the centerline. The third (self-similar) region begins at the distance from the die face at which the velocity maximum occurs on the centerline (see also Figure 7 in ref 12). The recirculation area in the first region of flow has a flow pattern that is the reverse of what is desired for melt blowing. The negative centerline velocity creates a drag force toward the die face, instead of a drag force in the direction of polymer flow. For a single annular jet, the recirculation area was significantly weakened by the presence of a fiber.26 Similarly, cases 4-6 had a much weaker recirculation area when compared to that present in the case 1 simulation. (These cases were compared because they all had the same Vjo.) Thus, the presence of a fiber is helpful in reducing negative centerline velocity. When a fiber is not present, the third region of the flow field is characterized by the maximum velocity occurring on the centerline. For a slot die, if a fiber is present, the maximum air velocity cannot occur along the centerline (although this maximum will still occur on the centerline between the fibers). The locations of the maximum velocity, xmax on the y ) 0 plane, for case 1 (air only) and for cases 4-6 (which all have the same nominal air velocity as case 1) are shown in Figure 6. These locations were determined by examining the velocity on lines beginning at the fiber edge (at y ) 0 and a constant z-position) to find the distance from x ) 0 where the maximum velocity occurred. Close to the die face, the location of the maximum velocity is governed by the location of the air jet. All simulations in Figure 6 show that the maxima move closer to the centerline, and all simulations are similar until z ≈ 2 mm. At z > 3 mm, the maximum for case 1 is located on the centerline, while the maxima for cases 4-6 are near x ) 0.3 mm. The location of the maximum for cases 4-6 appears to be constant and the same for all cases, indicating that the third, self-similar flow

940

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008

Figure 7. Dimensionless velocity for the centerline of case 1 and at y ) Lc/2 and y ) 0 for case 4.

Figure 8. Contour plots of z-velocity for case 4 at z ) 1, 2.5, 5, and 10 mm below the die face. The plot dimensions are to scale, and the diameter of the fiber at the top figure is Dc ) 0.475 mm.

region starts at z ≈ 3 mm. However, when the simulation results for the entire computational domain are closely examined, there is a slight spreading rate of xmax. At positions of z > 10 mm, a constant spreading of xmax for cases 2-7 can be described by the following equation (with both xmax and z in millimeters):

xmax ) S′z

(5)

For all cases, S′ ≈ 0.005. Also, for all cases, xmax < 1 mm throughout the entire 50 mm computational domain. The first goal of this paper is to investigate how the presence of a fiber affects the air flow field. Since the fiber moves along the centerline (x ) 0) in melt blowing, this is an area of interest. Figure 7 shows the dimensionless centerline velocity for case 1 (no fiber) and at two different y-locations for case 4. The first centerline location for case 4 is halfway between two fibers (represented by the periodic boundary in the simulations). Since the stress on the fiber is determined by the gradients close to the fiber, it is also of interest to look at the maximum velocity close to the fiber. The velocity at the points with y ) 0 (center of the domain) and x ) xmax is plotted in Figure 7. Although the position of xmax is difficult to measure experimentally, CFD allows fast and simple determination of xmax. As shown in Figure 6, the maximum velocity does not occur along the centerline at y ) 0 due to the presence of the fiber. The velocity profiles for case 1 and both positions for case 4 show significant differences. The case 4 profile at y ) Lc/2 has the highest maximum, which can be attributed to the fact that space available for the air flow is reduced due to the presence of the fibers. The case 4 profile at y ) 0 has the lowest maximum because the presence of a slow moving fiber, which moves slower than the air, leads to a drag force that reduces the air velocity. The x ) xmax position is very close to the fiber, and the difference between the air velocity and the fiber velocity results in the attenuating drag force. This velocity gradient is of interest because this is the air velocity that the fiber “feels”, and this gradient plays a critical role in the mechanism of the melt blowing process. Figure 7 shows that the velocity profiles of the air are nonuniform and are altered due to the presence of a fiber. Figure 8, which shows the z-velocity contours for case 6, illustrates the development of the flow field with the inclusion of a fiber. At z ) 1 mm the two bands representing the highest air speeds show that the velocity maxima are determined by the location of the jet (also see Figure 6); the two jets have not merged together. The second contour is at z ) 2.5 mm below

Figure 9. Dimensionless centerline velocity for all cases at y ) Lc/2. The R2 ) 0.9835 for the fitted equation shown on the figure.

the die face. The two jets have merged together as much as possible with the fiber between them; the maximum air velocity occurs on either side of the fiber at y ) 0. A strong attenuating force is present at this location (see Figure 14), but the force on the fiber is not uniform due the nonuniform velocity profiles. The next contour, at z ) 5 mm, shows that the air velocity maxima are now located halfway between the fibers, at y ) (Lc/2. This is the location of the velocity maxima for the computational domain beyond this point, as is also shown in the final contour plot at z ) 10 mm. Figure 8 shows the development of the velocity field when a fiber is included. The velocity gradients around the fiber are not uniform. Figure 9 shows the dimensionless centerline velocity profiles halfway between two fibers, at y ) Lc/2. The location and magnitude of the velocity maxima for cases 2-7 are similar, but are different from that of case 1 (air-only simulation). The air is being squeezed by the presence of the fiber close to the die face (where the diameter is the greatest), which leads to higher air velocity. The velocity for case 7 decays more rapidly than for the other cases. This is the case with the lowest air flow rate; the slow moving polymer is creating a slowing drag force on the air flow. Far from the die, the velocity decays of cases 2-6 are close to each other. Figure 10 shows the dimensionless velocity on the y ) 0 plane for all cases. For cases 2-7, this velocity occurs at x ) xmax. When compared to the air velocities close to the fiber, the air-only simulation

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 941

Figure 10. Dimensionless centerline velocity for all cases at y ) 0 and x ) xmax. The R2 ) 0.9835 for the fitted equation shown in the figure. Table 3. Velocity Decay Constants for y ) 0 and x ) xmax Starting at z > 10 mm case

a

b

R2

2 3 4 5 6 7

1.2613 1.201 1.201 1.0243 1.0566 0.9938

-0.5813 -0.5516 -0.5434 -0.509 -0.5007 -0.5405

0.9979 0.9989 0.9993 0.9975 0.9993 0.9985

Figure 11. Turbulence intensity as a percentage of Vjo on the centerline at y ) Lc/2.

reaches a higher maximum because the air is not being slowed by a fiber (which is moving slowly close to the die face). It has been shown6 that the velocity decay from a single planar jet can be described using the equation

Vo ) c1Vjoz-1/2

(6)

where Vo is the air velocity along the centerline, Vjo is the nominal velocity at the die face, and c1 is an empirical constant. Harpham and Shambaugh23 measured the centerline velocity decay for a blunt slot die and developed the following correlation:

Vo ) 1.4Vjo*z-0.610

(7)

Except for the slight difference in the exponent, eq 7 is the same as eq 6. Although the profiles for cases 2-7 in Figure 10 are not directly on the centerline, a power law correlation similar to eqs 6 and 7 can be used to describe the decay for these cases for z > 10 mm. Table 3 gives the values of c1 for all cases that included a fiber; the c1 values range from 0.9938 to 1.2613. In general, the value of this constant decreases as the momentum flux ratio increases. The exponent for cases 2-7 is also given in Table 3. In theory, this exponent for a single planar jet is -0.5 (eq 6). The value of the exponent for cases 2-7 ranges between -0.5813 and -0.5007. These values are between the theoretical value for a single plane jet and the experimentally measured decay (see eq 7) for a blunt melt blowing die when no fiber is present. If the best fit constants and exponents for cases 2-7 are averaged, the following equation results:

Vo ) 1.123Vjo*z-0.5378

(8)

The R2 value for eq 8 is 0.9835. Therefore, this equation can be used to accurately predict the velocity at y ) 0 and x ) xmax. This is the velocity of the air that is closest to the fiber,

Figure 12. Contour plots of turbulence intensity as a percentage of Vjo. These contour plots are for xy-planes located at z ) 1.2 mm and z ) 10 mm for cases 1 and 6. The plot dimensions are to scale, and the diameter of the fiber in the top figure is Dc ) 0.475 mm.

and it is the air velocity responsible for the attenuation of the fiber diameter. The trendline of eq 8 is included in Figure 10. 3.3. Fluctuating Velocity Field. Simulations of an annular jet with and without a fiber revealed that turbulence in the air flow field was decreased due to the presence of a fiber.26 For the slot jet in the present study, Figure 11 shows the turbulence intensity, q, along the periodic centerline, as a percentage of Vjo, for all simulations in Table 1. Case 1, which does not include a fiber, has the highest turbulence intensity maximum (although the location of this maximum is the same for all simulations). Throughout the entire computational domain, case 1 has higher turbulence intensity. Therefore, even at the centerline location farthest from the fibers, the turbulence intensity is dampened significantly. In addition, examination of cases 2-7 reveals that, in the far field (where the largest differences in turbulence intensity profile occur), a higher momentum flux ratio leads to lower turbulence intensity. Figure 12 shows the contour plots of turbulence intensity for cases 1 and 6 at the positions of z ) 1.2 mm and z ) 10 mm below the die face. The location of z ) 1.2 mm corresponds to the distance from the die face where all cases exhibit a maximum (see Figure 11). For case 1 at z ) 1.2 mm, the turbulence intensity is clearly the highest at the center of the domain. Turbulence intensity for case 6 is also at a maximum on the centerline halfway between the fibers, where there is the least effect from the dampening presence of the fiber. However, it is clear that, even at the location of maximum turbulence intensity, case 6 has smaller velocity fluctuations than case 1. At z ) 10

942

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008

Figure 13. Dimensionless turbulence dissipation rate for all cases on the centerline at y ) Lc/2.

mm, where the flow field is self-similar, the turbulence intensity remains the highest along the centerline for case 1. However, for case 6 the turbulence intensity is dampened and is not on the centerline. Figure 13 shows the dimensionless turbulence dissipation rate for cases 1-7 at y ) Lc/2. All of the cases exhibit local maxima near z/h ≈ 0.3 mm. However, case 1 has the lowest peak at this point. Because the rate of turbulence dissipation is higher for the cases that include a fiber, these cases have less turbulent flow fields. Similar to the results for the annular air jet, the velocity fluctuations are decreased due to the presence of a fiber, while the dissipation of turbulence is increased by the fiber.26 This is because the fiber provides a solid boundary in an area of the flow field where the fluctuations were high in the aironly case. This solid wall has a no-slip boundary, which leads to zero fluctuations on the solid and small fluctuations close to this boundary. For optimal melt blowing, the air would act only to create a drag force on the fiber. Turbulence fluctuations do not help to attenuate the fiber and, thus, are a waste of energy. 3.4. Stresses on the Fiber Edge. The wall shear stress on the fiber is a result of the difference between the air and fiber speeds. The wall shear stress can be defined as follows:

(

τw ) Fν

)

∂Va,eff,PAR ∂r

r)rf(z)

(9)

Figure 14 shows the wall shear stress on the surface of the fiber as predicted from the simulations. Since the velocity flow field surrounding the fibers is not uniform (see Figure 8), the wall shear stress at the same z-position, but at different radial locations of the fiber, is also different. Figure 14 shows the averaged wall shear stress at different distances below the die face. The average wall shear stress changes depending on both the initial air velocity and the polymer flow rate, since the stress is related to the difference in the velocity of the air and fiber. The wall shear stress from case 2 has the highest maximum, because case 2 has the largest difference between the air and fiber speed. Cases 4-6 all have the same air flow rates, so the difference in the shear stress for these cases can be explained by the different fiber speeds. Out of these three cases, case 6 has the largest polymer flow rate, which leads to the lowest fiber speed. Therefore, case 6 has a higher wall shear stress maximum, while case 4 (with the lowest polymer flow rate) shows a higher wall shear stress maximum compared to cases 5 and 6. Finally, case 7, with the lowest air flow rate, has the lowest wall shear stress maximum.

Figure 14. Averaged wall shear stress on the fiber edge for cases 2-7. Table 4. Wall Shear Stress Deviation case

maximum coefficient of variation

average coefficient of variation

2 3 4 5 6 7

0.562 0.572 0.520 0.542 0.711 0.517

0.117 0.098 0.119 0.116 0.102 0.087

As mentioned previously, the wall shear stress is not constant at each z-location. To quantify the variance in the wall shear stress at any given z-position, Table 4 gives the maximum coefficient of variation (the standard deviation at each z-position divided by the average wall shear stress at that position) and the z-position where this maximum occurs for each of the cases. To calculate the maximum coefficient of variation, the wall shear stress reported at six radial positions (corresponding to six different computational cells of the computational grid) was used. The maximum coefficient of variation for all cases is close to the die face at z < 0.5 mm. Table 4 also provides the average coefficient of variation for the wall shear stress. Matsui31 developed the following empirical relation to calculate the friction factor between the air and fiber in melt spinning:

Cf ) β(ReDP)-n

(10)

The Reynolds number based on fiber diameter, ReDP, is defined as

ReDP ) FaVa,eff,PARdf/µa

(11)

The exponent in eq 10, n, was determined by Matsui to be 0.61.31 Later, Majumdar and Shambaugh32 experimentally determined that, for the higher velocity conditions of melt blowing, β ) 0.78 was a better value to use. Using the friction coefficient as defined in eq 10, Marla and Shambaugh17 calculated the drag force on the fiber as follows:

FPAR ) Cf(1/2)FaVa,eff,PAR2πdf Lf

(12)

Va,eff,PAR is the difference between the air and fiber speeds in the z-direction. The experiments of Majumdar and Shambaugh,32 as well as the modeling work of Marla and Shambaugh,17 are based on the assumption that the air flow field is not changed by the presence of a fiber. Using CFD, this assumption can be

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 943 Table 5. Comparison of Calculated Drag Force and Simulation Results cases 1, 4 1, 5 1, 6

percent difference

b

5.47 18.66 34.25

0.825 0.956 1.186

Table 6. Spreading Rates for z g 30 mm case

Sy)0, x)xmax

Sy)Lc/2

1 2 3 4 5 6 7

N/A 0.0922 0.1096 0.1296 0.1295 0.1197 0.1255

0.1091 0.1357 0.1159 0.1398 0.1407 0.1449 0.1369

tested, since the CFD package reports the wall shear stress. The centerline velocity from case 1 (air only) was used with the fiber speed from cases 4-6 (that have the same Vjo as case 1) to calculate Va,eff,PAR. Then eqs 9-11 were used to determine the wall shear stress, which is the drag force divided by the fiber surface area. As mentioned above, the value of β ) 0.78 was determined experimentally. Table 5 gives β′, which is the value of β that must be used in eqs 9-11 so that the wall shear stresses are the same as those predicted from the CFD simulation. For the lowest momentum flux ratio (case 4), the β′ is only 5.47% higher than β ) 0.78. However, this percent difference rises to 34.25% at the highest flux ratio (case 5). This suggests that eqs 9-11 work better when the fiber flow rate is lower. 3.5. Jet Spreading Rate. The spreading rate for a jet can be described as follows:6

x1/2/h ) S(z/h) + A

(13)

where x1/2 is the distance from the centerline where the air velocity is equal to half the centerline velocity. The constant A is related to the virtual origin of the jet. For a slot melt blowing die with a blunt tip, Harpham and Shambaugh23 found that S ) 0.118. When a fiber is present at the center of the domain, then the location of xmax cannot occur at the centerline (xmax * 0). Then, the following equation can be used to calculate the spreading rate:26

(x1/2 - xmax) ) S(z/h) + B h

(14)

Table 6 gives the spreading rates for all cases at two positions. The first position is halfway between the fiber capillaries at y ) Lc/2, which is calculated using eq 13. The second position is at y ) 0, where the fiber is present. For this position, eq 14 must be used to determine the spreading rate. For all cases, the spreading rate at y ) Lc/2 is higher than at y ) 0. 3.6. Temperature Field. To more closely model industrial melt blowing conditions, all simulations discussed in this paper were completed at non-isothermal conditions. For all simulations, the initial temperature of the fiber was 295 °C, the initial excess temperature of the air (Θjo) was 309 °C, and the temperature of the air entrained by the jet through the pressure outlet boundaries was 21 °C. Figure 15 shows the excess temperature at y ) Lc/2 for all cases. Case 1 (air only) had the most rapid temperature decay. This is because the fiber present in cases 2-7 was hot and released heat to the air. The temperature depended on the momentum flux ratio. The

Figure 15. Temperature decay on the centerline at y ) Lc/2 for all cases.

Figure 16. Contours of excess temperature for case 6 at z ) 1 mm and z ) 2.5 mm. The plot dimensions are to scale, and the diameter of the fiber at the top figure is Dc ) 0.475 mm. Table 7. Excess Temperature Decay Constants for at z g 5 mm case

c

d

R2

1 2 3 4 5 6 7

0.9961 1.1053 1.0937 1.0085 1.0185 1.0196 0.9981

-0.5911 -0.4715 -0.4566 -0.3735 -0.3630 -0.3374 -0.2882

0.9831 0.9987 0.9959 0.9994 0.9992 0.9989 0.9956

simulation with the highest momentum flux ratio, case 7, also exhibited the highest temperature because the fiber was larger and thus provided more thermal energy to heat the air. Also, for case 7, the air was moving more slowly, and thus there was more time for energy to transfer from the fiber to the air. The mean velocity, fluctuations, and spreading rate have been discussed in the previous sections for the positions y ) 0 and y ) Lc/2. Figure 16 shows the contour plots of the excess temperature at z ) 1 mm and z ) 2.5 mm. At these positions, as well as throughout the computational domain, the excess temperature field does not change significantly between different y-positions. Therefore, correlations for temperature decay at y ) Lc/2 will also apply at y ) 0 and x ) xmax. Experimental23,24 and computational10 studies have shown that the temperature decay from a non-isothermal jet can be described by a relation of the form

Θ/Θjo ) c*(z/h)d

(15)

For cases 1-7 the excess temperature decay constants for z g 5 mm are given in Table 7. Harpham and Shambaugh23 experimentally determined c and d to be 1.2 and -0.615, respectively. As the momentum flux ratio increased, the decay exponent, d, decreased because the fiber released more heat to the air, which slowed the excess temperature decay.

944

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008

4. Conclusions Computational fluid dynamics were used to study the interactions between air and a fiber during the melt blowing process for an Exxon slot die. A rheological (polymer-side) model was used to predict the fiber velocity, temperature, and diameter. Then, these predictions were used as boundary conditions for the simulation of air flow around the fiber. The presence of the fiber was found to have an important impact on the air flow. The maximum of the centerline velocity of the air is slightly increased relative to the air-only case at y ) Lc/2, which corresponds to a position halfway between two fibers. At y ) 0, the decay of the maximum air velocity for fiberinclusive simulations can be modeled with a single curve. Since the simulation results showed the air profile to be nonuniform, the drag force depends on the radial position on the fiber edge. The turbulence intensity is dampened by the presence of a fiber, while the dissipation of turbulence is increased. Lower turbulence intensity means less velocity fluctuations, which allows for higher air velocities to be used. Thus, the presence of the fiber actually improves air conditions for the process of melt blowing. The jet spreading rate is increased by the presence of a fiber, and the spreading rate halfway between fibers is higher than the spreading rate at y ) 0. As the momentum flux ratio increases, the rate of the excess temperature decay of the air is decreased. Acknowledgment This work was supported by an NSF GOALI grant (DMII0245324). The support of 3M, Procter & Gamble, and ConocoPhillips is gratefully acknowledged, as well as the permission by Fluent Inc. to use the FLUENT 6.2 software with an educational license. Nomenclature A, B, C ) domain edges to indicate different cell spacing (see Figure 2) A ) constant in spreading rate equation (see eq 13) b ) jet slot width at die face, mm B ) constant in spreading rate equation (see eq 14) Cf ) drag coefficient cp ) air heat capacity, J/K C1, C2 ) parameters for the dissipation equation in the RSM model c1 ) empirical parameter for centerline air velocity decay (eq 6) c, d ) parameters in the excess temperature decay equation (eq 15) Dc ) diameter of polymer capillary, mm df ) fiber diameter, mm E ) energy, J Gk ) production of turbulent kinetic energy due to buoyancy h ) distance from the beginning of one air jet to the end of the next air jet, mm k ) turbulent kinetic energy (1/2uiui), m2/s2 kc ) thermal conductivity, J/(m‚s‚K) Lc ) domain width in y-direction, mm mf ) mass flow rate of polymer, g/min ma ) mass flow rate of air, g/min p ) pressure fluctuation, Pa Prt ) turbulent Prandtl number, Prt ) 0.85 qj ) momentum flux ratio (see eq 1) q ) turbulence intensity, 1/3(uiui)1/2/Vjo × 100%

Qair ) air flow rate, slpm r ) spatial coordinate, mm Re ) Reynolds number of the air flow based on the individual jet width ReDP ) Reynolds number based on Va,eff,PAR and fiber diameter (eq 11) rf ) fiber radius, µm or mm R2 ) coefficient of determination for the evaluation of regression analysis S ) momentum spreading rate (eq 14) S′ ) spreading rate of xmax (eq 5) t ) time, s T ) temperature, K Ui ) mean velocity in the ith direction, m/s ui ) velocity fluctuation in the ith direction, m/s uiuj (i,j) component of the Reynolds stress divided by the fluid density, m2/s2 V ) air velocity in the z-direction, m/s Va,eff,PAR ) difference between fiber and air z-velocity, m/s Vf ) fiber velocity in the z-direction, m/s Vf,die ) fiber velocity at the die face, m/s Vjo ) nominal discharge velocity defined as the volumetric air flow divided by the area available for flow, m/s Vo ) centerline air velocity, m/s xk ) kth Cartesian direction, mm x ) spatial coordinate, mm xmax ) the location of the maximum velocity on the y ) 0 plane, mm y ) spatial coordinate, mm z ) spatial coordinate, mm Greek Characters β ) empirical constant for calculating friction factor (eq 10) δkj ) Kroeneker’s delta  ) dissipation rate of turbulent kinetic energy, m2/s3 θ ) angle between air jet and die face, rad Θ ) excess air temperature defined as the difference between the ambient temperature and the air jet temperature, °C Θjo ) initial air excess temperature, °C µ ) viscosity, kg/(m‚s) µa ) air viscosity, kg/(m‚s) µt ) turbulent viscosity, kg/(m‚s) ν ) kinematic viscosity, m2/s F ) density, kg/m3 Fa ) air density, kg/m3 Fa,die ) air density at the die face, kg/m3 Ff,die ) polymer density at the die face, kg/m3 σ ) turbulent Prandtl number for the kinetic energy dissipation τw ) wall shear stress, Pa Literature Cited (1) Wente, V. A. Manufacture of Superfine Organic Fibers, Report PB111437, NRL-4364, April 15, 1954; Office of Technical Services, U.S. Department of Commerce: Washington, D.C., 1954. (2) Wente, V. A. Superfine Themoplastic Fibers. Ind. Eng. Chem. 1956, 48, 1342. (3) Buntin, R. R.; Keller, J. P.; Harding, J. W. Nonwoven Mats by Melt Blowing. U.S. Patent 3,849,241, Nov. 19, 1974. (4) Shambaugh, R. L. A Macroscopic View of the Melt Blowing Process for Producing Microfibers. Ind. Eng. Chem. Res. 1988, 27, 2363-2372. (5) Bradbury, L. J. S. The Structure of a Self-Preserving Turbulent Plane Jet. J. Fluid Mech. 1965, 23 (1), 31-64. (6) Pope, S. B. Turbulent Flows; Cambridge University Press: Cambridge, U.K., 2000; pp 134-135. (7) Heskestad, G. Hot-Wire Measurements in a Plane Turbulent Jet. J. Appl. Mech. 1965, 32, 721-734.

Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 945 (8) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Analysis of a Melt-Blowing Die: Comparison of CFD and Experiment. Ind. Eng. Chem. Res. 2002, 41, 5125-5138. (9) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Effects of Die Geometry on the Flow Field of the Melt Blowing Process. Ind. Eng. Chem. Res. 2003, 42, 5541-5553. (10) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Effects of Temperature and Geometry on the Flow Field of the Melt Blowing Process. Ind. Eng. Chem. Res. 2004, 43, 4199-4210. (11) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Using Computational Fluid Dynamics to Simulate Flow Fields from Various Melt Blowing Dies. Int. NonwoVens J. 2005, 14 (1), 2-8. (12) Moore, E. M.; Shambaugh, R. L.; Papavassiliou D. V. Analysis of Isothermal Annular Jets: Comparison of CFD and Experimental Data. J. Appl. Polym. Sci. 2004, 94 (3), 909-922. (13) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Analysis of the Temperature Field from Multiple Jets in the Schwarz Melt Blowing Die Using Computational Fluid Dynamics. Ind. Eng. Chem. Res. 2006, 45, 5098-5109. (14) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Analysis of Multiple Jets in the Schwarz Melt-Blowing Die Using Computational Fluid Dynamics. Ind. Eng. Chem. Res. 2005, 44, 8922-8932. (15) Uyttendaele, M. A. J.; Shambaugh, R. L. Melt Blowing: General Equation Development and Experimental Verification. AIChE J. 1990, 36 (2), 175-186. (16) Rao, R. S.; Shambaugh, R. L. Vibration and Stability in the Melt Blowing Process. Ind. Eng. Chem. Res. 1993, 32, 3100-3111. (17) Marla, V. T.; Shambaugh, R. L. Three-Dimensional Model of the Melt-Blowing Process. Ind. Eng. Chem. Res. 2003, 42, 6993-7005. (18) Marla, V.; Shambaugh, R. L.; Papavassiliou, D. V. Modeling the Melt Blowing of Hollow Fibers. Ind. Eng. Chem. Res. 2006, 45, 407-415. (19) Yin, H.; Yan, Z.; Ko, W.-C.; Bresee, R.R. Fundamental Description of the Meltblowing Process. Int. NonwoVens J. 2000, 9 (4), 25-28. (20) Wang, K.; Ke, Q. Computational Simulation of the Fiber Movement in the Melt-Blowing Process. Ind. Eng. Chem. Res. 2005, 44, 3912-3917.

(21) Wang, K.; Ke, Q. Empirical Formulas for Distributions of Air Velocity/Temperature Along the Spinline of a Dual Slot Die. Polym. Eng. Sci. 2005, 45 (8), 1092-1097. (22) Hietel, D.; Marheineke, N. Modeling and Numerical Simulation of Fiber Dynamics. Proc. Appl. Math. Mech. 2005, 5, 667-670. (23) Harpham, A. S.; Shambaugh, R. L. Flow Field of Practical Dual Rectangular Jets. Ind. Eng. Chem. Res. 1996, 35, 3776-3781. (24) Harpham, A. S.; Shambaugh, R. L. Velocity and Temperature Fields of Dual Rectangular Jets. Ind. Eng. Chem. Res. 1997, 36, 3937-3943. (25) Tate, B. D.; Shambaugh, R. L. Temperature Fields below Melt Blowing Dies of Various Geometries. Ind. Eng. Chem. Res. 2004, 43, 54055410. (26) Krutka, H. M.; Shambaugh, R. L.; Papavassiliou, D. V. Effects of the Polymer Fiber on the Flow Field from an Annular Melt Blowing Die. Ind. Eng. Chem. Res. 2007, 46, 655-666. (27) Fluent User’s Guide; Fluent Inc.: Lebanon, NH, 2007. (28) Gambit Modeling Guide; Fluent Inc.: Lebanon, NH, 2007. (29) Schetz, J. A.; Padhye, A. Penetration and Breakup of Liquids in Subsonic Airstreams. AIAA J. 1977, 15 (10), 1385-1390. (30) Nasr, A.; Lai, J. C. S. Comparison of Flow Characteristics in the Near Field of Two Parallel Plane Jets and an Offset Plane Jet. Phys. Fluids 1997, 9 (10), 2919-2931. (31) Matsui, M. Air Drag on a Continuous Filament in Melt Spinning. Trans. Soc. Rheol. 1976, 20, 465-473. (32) Majumdar, B.; Shambaugh, R. L. Air Drag on Filaments in the Melt Blowing Process. J. Rheol. 1990, 34 (4), 591-601.

ReceiVed for reView June 25, 2007 ReVised manuscript receiVed October 4, 2007 Accepted October 10, 2007 IE070871I