# Optimization-based design of polymer sheeting dies using generalized Newtonian fluid models.

INTRODUCTIONThe analysis and design of the polymer extrusion process used to produce thin sheets and films have received much attention over the last few decades. While significant advances have been made, complex constitutive relations and demanding performance criteria continue to present a challenge to the successful design and operation of a sheeting die. Of particular interest is the design of a die cavity geometry having a minimal inlet pressure while delivering a thin layer of polymer melt with a uniform velocity over the entire width of the die exit.

Closed-form sheeting die design equations have been derived assuming one-dimensional flow of isothermal Newtonian and power-law fluids in the manifold and land regions of a sheeting die [1-5]. These early developments have been applied to T-dies, fish-tail dies, and coat-hanger

dies in which the cross-section geometry of the distribution manifold is defined over the width of the die to yield a uniform exit velocity. Sheeting dies have also been designed based on one-dimensional flow for slowly reacting polymeric liquids [6]. More recently, Liu et al. [7] developed a similar but more general approach that is capable of incorporating Newtonian, power-law, modified Bingham, and Ellis fluid models. The resulting closed form die design equations provide a computationally efficient approach for defining sheeting die cavity geometries, but assume a constant wall shear rate and are limited to a finite set of geometric configurations.

The flow of polymer melt in sheeting dies has also been simulated using three-dimensional flow equations and numerical approaches such as the finite element method to solve for the pressures, velocities, and temperatures in the polymer melt [8-10]. Various generalized Newtonian fluid (GNF) models (see, e.g., Tadmor and Gogos [11] or Bird et al. [12]) have been included in these simulations, which have been shown to compare well with experimental results by Dooley [8]. Simplifying assumptions such as constant temperature and unidirectional fully developed flow, which serve as a basis for the one-dimensional techniques described above, are avoided. Instead, a complete description of the melt flow is obtained for arbitrary die cavity geometries that are not restricted to simple T, fish-tail, or coat-hanger shapes. Recently, numerical optimization methods have been combined with three-dimensional melt flow simulations, which provides a powerful technique for designing die cavities [13, 14]. This type of design methodology allows for a more general definition of the die cavity and may include other important performance metrics such as die pressure drop. Unfortunately, the computational effort required to obtain the more accurate three-dimensional solutions can be excessive when used in an iterative design procedure, particularly when applied to detailed industrial designs.

Two-dimensional (often referred to as 2.5 D) simulations provide a computationally efficient alternative to three-dimensional methods, while avoiding many of the simplifying assumptions and geometric restrictions of existing one-dimensional design equations (see, e.g., Refs. 15-20). Sartor [15] was perhaps the first to combine numerical optimization with a two-dimensional flow network analysis to iteratively solve the die design problem using a power-law fluid model. Smith et al. [16, 17] developed a finite element simulation based on the Hele-Shaw flow approximation [21] similar to that used in injection molding simulations [22] to compute optimal die cavity geometries. Unfortunately, these designs were also limited to isothermal Newtonian and power-law fluids. Constraints were defined to measure the die's exit velocity and design sensitivities were developed in Refs. 16 and 17, which greatly reduced the computational effort required in the design calculations. This approach was extended to dies having a uniform exit velocity and residence time; examples include coat-hanger dies [17] and dies with a general die cavity height distribution [19]. More recently, this die design approach was extended to address variability during operation [23]. Die cavity deformations were incorporated in the die design method by Bates et al. [24], which computed optimal restrictor geometries using power-law fluids and the Hele-Shaw flow approximation.

This article develops and applies an optimization-based approach that computes die cavity geometries to deliver a uniform exit velocity at a specified flow rate and minimum die pressure drop. Isothermal flow simulations are performed with the finite element method, in which melt flow calculations are simplified with the Hele-Shaw flow approximation. Of particular interest here are the analysis and design sensitivity analysis for polymer melt flows represented by GNF models. Integrals through the die cavity thickness for flow conductance are evaluated numerically to accommodate a Bingham, Ellis, Cross, and Carreau-Yasuda viscosity model. The system of nonlinear melt flow equations is solved via a Newton-Raphson iteration that facilitates efficient computation of the design sensitivities. Similarities between tangent matrix formulations and design sensitivities are exposed. The proposed design method provides a computationally efficient means to compute optimal die cavity shapes for various purely viscous fluid models. Examples are given that compare designs computed when different polymer melt models are employed. Pareto optimal designs (see, e.g., Ref. 25) are also computed to expose the pressure drop versus flow rate relationship for the given die design parameterization.

SHEETING DIE DESIGN

In the die design problem considered in this paper, the pressure drop across the die and the die exit velocity variation define the success of a given die design. They are chosen since pressure drop determines the extruder size and power requirements, and die exit velocity variation influences the sheet thickness uniformity [17]. The goal of minimizing the inlet pressure is realized by varying the thickness distribution in the die cavity, while placing constraints on the die exit velocity variation and the flow rate of the polymer melt, and limiting the slope of die cavity surface. The die design optimization problem may be stated mathematically as

[FIGURE 1 OMITTED]

[min.[[phi][member of][R.sup.N]]] f([phi]) = [P.sub.in]

Such that [g.sub.1]([phi]) = [1/L][[integral].sub.l.sub.exit] ([[[bar.v]([phi])]/[[v.sub.a]([phi])]] - 1)[.sup.2] dx [less than or equal to] [[epsilon].sub.1]. (1)

[g.sub.2]([phi]) = ([[[v.sub.a]([phi])]/[bar.v.sub.p]] - 1)[.sup.2] [less than or equal to] [[epsilon].sub.2]

[g.sub.3.sup.K]([phi]) = [||[nabla][h.sup.K]([phi])||/||[nabla][h.sub.p]||] - 1 [less than or equal to] 0

The nonlinear constrained optimization (see, e.g., Arora [26]) given here is to minimize the inlet pressures [P.sub.in]. Design parameters are denoted by [phi] = {[P.sub.in], [h.sup.I]}, I = 1, 2 ... N - 1, where [h.sup.I] are the N - 1 nodal half-height design variables and N is the total number of design variables.

The constraint function [g.sub.1] in Eq. 1 measures the exit velocity variation and is included to obtain a uniform exit velocity to within the tolerance [[epsilon].sub.1]. This constraint is imposed to maintain the gap-wise average exit velocity [bar.v]([phi]) normal to the die exit in Fig. 1 at a value that approaches the average exit velocity [v.sub.a]([phi]) defined as

[v.sub.a]([phi]) = [1/L][[integral].sub.l.sub.exit] [bar.v]([phi]) dx (2)

where [l.sub.exit] denotes the die exit edge, and L is the total length along the die exit. A zero value of [g.sub.1] indicates the computed [bar.v]([phi]) equals [v.sub.a]([phi]) over the entire width of the die exit. The constraint function [g.sub.2] is imposed to set the flow rate through the die. In our example, the die exit has a uniform half-height [h.sub.exit] that does not change as a function of design. Therefore, when the average exit velocity [v.sub.a] equals the prescribed gap-wise average exit velocity [bar.v.sub.p], the die is operating at the desired total flow rate Q given as

Q = 2[h.sub.exit][[integral].sub.l.sub.exit] [bar.v.sub.p] dx. (3)

In the optimization problem in Eq. 1 the inlet pressure [P.sub.in] may vary with the design to achieve a desired flow rate [23]. This formulation provides a means to independently enforce exit velocity uniformity with [g.sub.1] and total flow rate with [g.sub.2], thus allowing more control over the design than in previous optimizations [17]. In the die cavity, the half-height h(x, y) may be arbitrary. Therefore, constraints [g.sub.3.sup.K], K = 1, 2 ... [N.sub.s], are imposed to restrict the slope of h(x, y) to within the prescribed value ||[nabla][h.sub.p]|| to avoid extensional flows that would otherwise invalidate the flow model assumptions.

It is important to note that the objective function f([phi]) and constraint [g.sub.3.sup.K]([phi]) are explicitly dependent on the design parameters, while the constraints [g.sub.1]([phi]) and [g.sub.2]([phi]) are implicitly dependent on the design parameters through the melt flow governing equations. This implicit dependence complicates the design sensitivity analysis. In the following sections, we present an analytical approach to sensitivity analysis and its implementation with the finite element method that may be used with various GNF models. Note that as with previous die design methods discussed above, our approach does not include the solution of the energy equation when computing new die geometries. This represents a simplification, rather than a limitation of the design methodology presented below.

MODELING AND SIMULATION

To optimize the die cavity geometry of a sheeting die such as that appearing in Fig. 1, a die flow simulation is required. The Hele-Shaw flow model [21] is developed from the principles of conservation of mass, momentum, and energy, where assumptions are made to reduce computational time and data requirements. The polymer melt is modeled as a GNF, and therefore is assumed to be inelastic and incompressible. In addition, inertial, body, and surface tension forces in the fluid are assumed to be negligible. Moreover, the pressure does not vary significantly in the direction normal to the plane of flow, and the velocity in the direction normal to the plane of flow is negligible compared with the in-plane velocities (see, e.g., Fig. 1). The die cavity thickness is assumed to be small compared to its in-plane dimensions and has little in-plane variation and all flow conditions are assumed to be symmetric with respect to the cavity mid-plane. This model is widely employed in injection and compression molding [21], and has also been applied to sheet extrusion dies [15, 17, 19, 20, 23, 24].

Governing Equations

Based on these assumptions, the mass and momentum conservation equations reduce to a single differential equation [21]

[nabla] * S[nabla]P = 0 (4)

where P is the pressure field over the flow cavity domain [OMEGA] in two-dimensions, and [nabla] is the gradient operator in the x - y plane. The boundary conditions in the plane of flow are P = [P.sup.p] on [partial derivative][[OMEGA].sup.p] and S([nabla]P * n) = [q.sup.p] on [partial derivative][[OMEGA].sup.q], where [q.sup.p] is the prescribed flow rate. Hence, the boundary of the flow domain [partial derivative][OMEGA] is divided into two complimentary subdomains: [partial derivative][[OMEGA].sup.p] and [partial derivative][[OMEGA].sup.q]. S is the flow conductance defined as an integral through the cavity thickness as [21]

S = [[integral].sub.0.sup.h] [[z.sup.2]/[[eta](z)]] dz. (5)

A residual R for the boundary value problem described in Eq. 4 is obtained via the method of weighted residuals (see e.g., Ref. 27) in the usual manner as [19]

R(P) = [[integral].sub.[OMEGA]] [nabla]w * S(P)[nabla]Pda - [[integral].sub.[partial derivative][OMEGA].sup.q] w[q.sup.p] ds (6)

where w is an arbitrary weighting function. The pressure field P that satisfies Eq. 6 is the solution to Eq. 4 if w vanishes on [partial derivative][[OMEGA].sup.p], and P satisfies the boundary conditions on [partial derivative][[OMEGA].sup.p] [17]. When non-Newtonian fluids are considered, S in Eq. 5 becomes a function of the pressure field and Eq. 6 is nonlinear in P, requiring iterative methods to compute a solution. In this research, the nonlinear residual Eq. 6 is solved via the Newton-Raphson method since it exhibits terminal quadrature convergence and is conducive to the design sensitivity analysis to follow. The tangent operator [[partial derivative]R]/[[partial derivative]P] operating on the increment [[DELTA]P] is obtained by differentiating Eq. 6 with respect to P as

[[[partial derivative]R(P)]/[[partial derivative]P]][[DELTA]P] = [[integral].sub.[OMEGA]][nabla]w * [S(P)[nabla][[DELTA]P] + [[[partial derivative]S(P)]/[[partial derivative]P]][[DELTA]P][nabla]P]da (7)

where we have assumed that [q.sup.p] in Eq. 6 is not a function of the pressure P. Note that Eqs. 6 and 7 may include any GNF model in which the material model enters the computations through S and [[[partial derivative]S]/[[partial derivative]P]].

Model Solution

The isoparametric finite element method (see, e.g., Ref. 27) is used to discretize the above equations when the two-dimensional flow domain [OMEGA] in Eqs. 6 and 7 is divided into multiple finite element domains. Note that the die cavity half-height h enters the calculations through the integral in Eq. 5. In the finite element analysis, the melt pressure field P and residual R in Eq. 6 become the nodal pressure vector P and residual vector R, respectively (see, e.g., Refs. 17 and 19). In a similar manner, the tangent operator [[partial derivative]R]/[[partial derivative]P] in Eq. 7 becomes the tangent matrix [[partial derivative]R]/[[partial derivative]P]. Once the residual R and tangent matrix [[partial derivative]R]/[[partial derivative]P] are computed, the nodal pressures may be evaluated using Newton-Raphson iteration as

[[[partial derivative]R([P.sub.i])]/[[partial derivative]P]] [[DELTA]P] = -R([P.sub.i]) (8)

where nodal pressures are updated as [P.sub.i + 1] = [P.sub.i] + [DELTA]P. Iterations are repeated until convergence is reached. Once the pressure solution is computed, the gap-wise average velocity vector [bar.v] is evaluated from

[bar.v] = - [S/h][nabla]P. (9)

Additionally, the scalar measure [dot.[gamma]] = [square root of ([1/2][[PI].sub.[dot.[gamma]]])] of the strain rate tensor [dot.[gamma]] (recall that [[PI].sub.dot.[gamma]] is the second invariant of the strain rate tensor [dot.[gamma]] as described in Ref. 11) for the Hele-Shaw flow approximation is computed at a distance z from the die cavity mid-plane as

[dot.[gamma]] = [z/[eta]]||[nabla]P|| (10)

where ||[nabla]P|| is the magnitude of the pressure gradient in the plane of flow.

Generalized Newtonian Fluids

The GNF model is widely used to represent the purely viscous non-Newtonian behavior of polymer melt flow (see, e.g., Refs. 11 and 12), where the shear stress tensor [tau] is proportional to the strain-rate tensor [dot.[gamma]], as in

[tau] = [eta]([dot.[gamma]])[dot.[gamma]]. (11)

Expressions for several GNF models are given in Table 1, where model parameters are described in the literature (see, e.g., Refs. 11 and 22). When an isothermal power-law fluid is employed, the viscosity [eta] = m[dot.[gamma].sup.n-1] from Table 1 may be substituted into Eq. 10 to obtain an explicit expression for [eta] as a function of z. This result may be combined with Eq. 5 to evaluate the flow conductance S as [17]

[TABLE 1 OMITTED]

S(P) = [[h.sup.(1/n)+2]/[[m.sup.(1/n)]([1/n] + 2)]]||[nabla]P||[.sup.(1/n)-1]. (12)

Consequently, the term [[partial derivative]S]/[[partial derivative]P] acting on the increment [[DELTA]P] in the tangent operator [[partial derivative]R]/[[partial derivative]P] in Eq. 7 is obtained by differentiating Eq. 12 as [17]

[[[partial derivative]S]/[[partial derivative]P]][[DELTA]P] = [[S([1/n] - 1)]/[||[nabla]P||[.sup.2]]][nabla]P * [nabla][[DELTA]P] (13)

where we note that an analytical expression for S(P) and [[partial derivative]S]/[[partial derivative]P] are possible since m and n are independent of z for homogeneous isothermal flow. Note that Eqs. 12 and 13 reduce to S = [h.sup.3]/[3[mu]] and [[partial derivative]S]/[[partial derivative]P] = 0, respectively, for a Newtonian fluid with viscosity [mu].

Since the power-law fluid model does not capture the near-constant viscosity at low strain rates, other more realistic GNF models such as the Carreau-Yasuda, Ellis, and Cross models are often used. Unfortunately, the simplicity of the Newtonian and power-law implementation described above is lost when more complex GNF models are employed since the viscosity [eta] cannot be written as a function of ||[nabla]P|| in the generalized Hele-Shaw analysis. As a result, an analytical expression for the flow conductance S in Eq. 5 does not exist. Therefore, we employ Gauss quadrature (see, e.g., Ref. 27) to numerically integrate Eq. 5 as

S = [NG.summation over (i=1)][c.sub.i][[z.sub.i.sup.2]/[[eta]([z.sub.i])]] (14)

where NG is the number of integration points (which is NG = 8 in our analyses), [c.sub.i] are Gauss weighting factors, and the viscosity [eta] is evaluated at selected [z.sub.i] locations using equation Eq. 10 and the selected GNF viscosity formula from Table 1.

[TABLE 2 OMITTED]

Similarly, an analytical expression for the derivative [[partial derivative]S]/[[partial derivative]P] in Eq. 7, such as that shown in Eq. 13 for the power-law fluid, does not exist for the Ellis, Cross, Bingham, and Carreau-Yasuda fluid models. Since the domain of integration in Eq. 5 does not depend on P, we may differentiate its integrand to give

[[[partial derivative]S]/[[partial derivative]P]][[DELTA]P] = -[[integral].sub.0.sup.h][[z.sup.2]/[[eta].sup.2]][[[partial derivative][eta]]/[[partial derivative]P]][[DELTA]P]dz. (15)

The viscosity derivative [[partial derivative][eta]]/[[partial derivative]P] is evaluated by differentiating Eq. 10 and the selected viscosity formula in Table 1 with respect to the pressure P, which yields

[[[partial derivative][eta]]/[[partial derivative]P]][[DELTA]P] = [beta][[[nabla]P * [nabla][[DELTA]P]]/[||[nabla]P||[.sup.2]]] (16)

where the viscosity derivative factor [beta] is defined for each of the GNF models considered in this study in Table 2.

To numerically evaluate the integrals in Eqs. 5 and 15, the viscosity [eta] must be evaluated at selected z locations in a manner similar to that shown in Eq. 14. Unfortunately, when Eq. 10 is substituted into any of the viscosity equations given in Table 1, we obtain a nonlinear equation in [eta], which must be solved numerically via an iterative procedure such as the Newton-Raphson method. To illustrate the resulting viscosity calculation procedure, we consider the Carreau-Yasuda viscosity model given in Table 1. To this end, a residual expression r([eta]) based on the viscosity [eta] as well as its first derivative r'([eta]) = [dr([eta])]/[d[eta]] are defined, respectively, as

r([eta]) = [[[eta] - [[eta].sub.[infinity]]]/[[eta].sub.0]] - (1 - [[[eta].sub.[infinity]]/[[eta].sub.0]]) [1 + ([lambda][dot.[gamma]])[.sup.a]][.sup.(n-1)/a] (17)

r'([eta]) = [1/[[eta].sub.0]] + [[(n - 1)([eta] - [[eta].sub.[infinity]])([lambda][dot.[gamma]])[.sup.a]]/[[eta][[eta].sub.0](1 + ([lambda][dot.[gamma]])[.sup.a])]] (18)

where [dot.[gamma]] is computed for a given [nabla]P at each height z from Eq. 10. We desire that r([eta]) = 0, so that at iteration i, the viscosity is updated via the Newton-Raphson method as [[eta].sub.i+1] = [[eta].sub.i] - [[r([[eta].sub.i])]/[r'([[eta].sub.i])]], where iterations continue until convergence is achieved.

Once the viscosity [eta] is computed at all of the specified z locations, the flow conductance S and its derivative [[partial derivative]S]/[[partial derivative]P] may be obtained via numerical integration of Eqs. 5 and 15, respectively. The residual R and tangent operator [[partial derivative]R]/[[partial derivative]P] can then be evaluated from Eqs. 6 and 7, respectively.

DESIGN SENSITIVITY ANALYSIS (DSA)

Gradient-based optimization algorithms [26] use design sensitivities to efficiently compute improved designs. Design sensitivities quantify the relationships between the design variables such as the die cavity half-height or inlet pressure parameters in [phi] and system performance measures such as f, [g.sub.1], [g.sub.2], and [g.sub.3.sup.K] in Eq. I. In this study, we derive expressions for the design sensitivities, which are evaluated with the finite element method in a manner similar to that described above for the flow solution. This analytical approach avoids the well-known inaccuracies and inefficiencies associated with commonly used finite difference procedures.

Sheeting Die DSA

Since the objective function f([phi]) and constraint [g.sub.3.sup.K]([phi]) in Eq. 1 are explicitly defined on the design [phi], the calculation of design sensitivities is straightforward. For the objective function f([phi]) in Eq. 1, the design sensitivities may be computed from

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (19)

Following the same procedure, we obtain the sensitivity for the constraint [g.sub.3.sup.K]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

where the half-height gradients [nabla]h are only evaluated at the K-th element of interest.

In this study, the constraint functions [g.sub.1]([phi]) and [g.sub.2]([phi]) in Eq. 1 are not defined explicitly in terms of the design variables in [phi] given above. Instead, these functions are implicitly defined through the pressure solution computed by solving Eq. 6 via the finite element method. To evaluate the implicit sensitivities, it is convenient to consider a general performance measure

F([phi]) = G(P([phi]), [phi]) (21)

where F, which may represent either one of the implicit constraints, is defined through a function G with implicit dependence on [phi] through the nodal pressure vector P and explicit dependence on the design [phi]. Assuming sufficient smoothness, the design sensitivity of F with respect to the design variable [[phi].sub.i] may be calculated via the adjoint variable method (see, e.g., Ref. 17) as

[DF]/[D[[phi].sub.i]] = [[[partial derivative]G]/[[partial derivative][[phi].sub.i]]] - [lambda] * [[[partial derivative]R]/[[partial derivative][[phi].sub.i]]] (22)

where the adjoint variable vector [lambda] is computed from the linear system of equations [17]

[[[partial derivative]R]/[[partial derivative]P]][.sup.T][lambda] = ([[partial derivative]G]/[[partial derivative]P])[.sup.T]. (23)

In the above, [[partial derivative]G]/[[partial derivative]P] is the adjoint load associated with the performance measure of interest. The adjoint variable method is computationally efficient since it only requires the evaluation of one linear adjoint problem defined by Eq. 23 for each implicit objective and constraint function. Furthermore, it employs the transpose of the tangent matrix [[[partial derivative]R]/[[partial derivative]P]][.sup.T] from Eq. 8 used in the Newton-Raphson iteration for the pressure solution, which significantly reduces the computational effort.

Since the constraints [g.sub.1] and [g.sub.2] in Eq. 1 are not explicitly dependent on the design variables considered in this study, for all [[partial derivative]G]/[[partial derivative][[phi].sub.i]] = 0, [[phi].sub.i], i = 1 ... N. However, these constraints do depend on [phi] through the pressure solution. Therefore, we differentiate [g.sub.1] and [g.sub.2] with respect to the pressure vector P to, respectively, obtain the right-hand-side of Eq. 23 as

[[partial derivative][g.sub.1]]/[[partial derivative]P] = [[integral].sub.l.sub.exit]2([[bar.v]/[v.sub.a]] - 1)[[[[partial derivative][bar.v]]/[[partial derivative]P]][1/[v.sub.a]] - [[bar.v]/[v.sub.a.sup.2]][[[partial derivative][v.sub.a]]/[[partial derivative]P]]]dx (24)

[[partial derivative][g.sub.2]]/[[partial derivative]P] = 2([[v.sub.a]/[bar.v.sub.p]] - 1)([1/[bar.v.sub.p]][[[partial derivative][v.sub.a]]/[[partial derivative]P]]). (25)

It should be noted that the integral of Eq. 24 is defined over the entire length along the die exit. The derivatives [[partial derivative][bar.v]]/[[partial derivative]P] and [[partial derivative][v.sub.a]]/[[partial derivative]P] are, respectively, computed by differentiating Eq. 9 as

[[[partial derivative][bar.v]]/[[partial derivative]P]][[DELTA]P] = [1/h] [[[[partial derivative]S]/[[partial derivative]P]][[DELTA]P]||[nabla]P|| + [S/||[nabla]P||][nabla]P * [nabla][[DELTA]P]] (26)

and Eq. 2 as

[[[partial derivative][v.sub.a]]/[[partial derivative]P]][[DELTA]P] = [1/L][[integral].sub.l.sub.exit][[[partial derivative][bar.v]]/[[partial derivative]P]][[DELTA]P] dx. (27)

Following finite element discretization, we obtain [[partial derivative][bar.v]]/[[partial derivative]P] and [[partial derivative][v.sub.a]]/[[partial derivative]P] from Eqs. 26 and 27, respectively. By substituting [[partial derivative][bar.v]/[partial derivative]P] and [[partial derivative][v.sub.a]]/[[partial derivative]P] in Eqs. 24 and 25, the derivatives [[partial derivative][g.sub.1]]/[[partial derivative]P] and [[partial derivative][g.sub.2]]/[[partial derivative]P] are obtained.

To compute the design derivatives [[partial derivative]R]/[[partial derivative][[phi].sub.i]] in Eq. 22, we consider the pressure and height design variables separately. In the case of [[phi].sub.i] = [P.sub.in], the design derivative [[partial derivative]R]/[[partial derivative][[phi].sub.i]] is computed using the tangent matrix [[partial derivative]R]/[[partial derivative]P] as

[[partial derivative]R]/[[partial derivative][P.sub.in]] = [[[partial derivative]R]/[[partial derivative]P]][[[partial derivative]P]/[[partial derivative][P.sub.in]]] (28)

where [[partial derivative]P]/[[partial derivative][P.sub.in]] is a zero vector with unity in the components associated with [P.sub.in]. Therefore, only those nodes associated with inlet pressure [P.sub.in] need to be considered when computing Eq. 28. Note that [[partial derivative]R]/[[partial derivative][P.sub.in]] is identical for any GNF model.

When [[phi].sub.i] defines a cavity half-height, Eq. 6 is differentiated with respect to [[phi].sub.i] as

[[partial derivative]R]/[[partial derivative][[phi].sub.i]] = [[integral].sub.[OMEGA]][nabla]w * [[[partial derivative]S]/[[partial derivative][[phi].sub.i]]][nabla]P da (29)

which may be evaluated with the finite element method to obtain an expression for [[partial derivative]R]/[[partial derivative][[phi].sub.i]]. Here we assume [q.sup.p] in Eq. 6 is not a function of the half-height parameter. Note that the evaluation of [[partial derivative]R]/[[partial derivative][[phi].sub.i]] in Eq. 29 applies to any viscosity model described in Table 1, as described below.

Once the terms [[partial derivative][g.sub.1]]/[[partial derivative]P], [[partial derivative][g.sub.2]]/[[partial derivative]P], and [[partial derivative]R]/[[partial derivative][[phi].sub.i]] are computed, the adjoint variable vector [lambda] is obtained with Eq. 23 and the design sensitivity [DF]/[D[[phi].sub.i]] follows from Eq. 22. In these computations, we must evaluate [[partial derivative]S]/[[partial derivative]P] and [[partial derivative]S]/[[partial derivative][[phi].sub.i]], which depend on the fluid viscosity model. The finite element expression [[partial derivative]S]/[[partial derivative]P] can be derived from Eq. 15. In the ensuing section, we evaluate [[partial derivative]S]/[[partial derivative][[phi].sub.i]] for various GNF models, i.e., power-law, Carreau-Yasuda, Ellis, Cross, and Bingham fluid.

DSA for Generalized Newtonian Fluids

For an isothermal power-law fluid, the flow conductance S in Eq. 12 is an explicit function of the half-height h([phi]) so that the design derivative [[partial derivative]S]/[[partial derivative][[phi].sub.i]] in Eq. 29 may be computed from

[[partial derivative]S]/[[partial derivative][[phi].sub.i]] = [[([1/n] + 2)S]/h][[[partial derivative]h]/[[partial derivative][[phi].sub.i]]] (30)

where S from Eq. 12 is employed to simplify the final expression. Note that Eq. 30 reduces to [[partial derivative]S]/[[partial derivative][[phi].sub.i]] = [[3S]/h][[[partial derivative]h]/[[partial derivative][[phi].sub.i]]] for a Newtonian fluid.

In the examples that follow, we define the cavity half-height h(x, y) over each element as h(x, y) = [n.Nh] where N contains the element interpolation functions and [n.h] is the nodal half-height vector. When a nodal half-height [h.sup.J] is defined by the design variable [[phi].sub.i], the height sensitivity [[partial derivative]h]/[[partial derivative][[phi].sub.i]] with respect to that nodal half-height design variable at node J, i.e., [[phi].sub.i] = [h.sup.J], is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

where [N.sup.I] is the I-th component of the interpolation function matrix N and [[delta].sub.IJ] is the Kronecker delta ([[delta].sub.IJ] = 1 for I = J and [[delta].sub.IJ] = 0 for I [not equal to] J). Here we only consider [n.sub.e] elements that contain the J-th node.

Unfortunately, when any of the other GNF models other than Newtonian or power-law is employed, the derivative [[partial derivative]S]/[[partial derivative][[phi].sub.i]] cannot be obtained analytically since an explicit expression for S does not exist. Therefore, the flow conductance sensitivity must be calculated numerically. The flow conductance S of Eq. 5 may be rewritten to express its dependence on [phi] as

S([phi]) = [[integral].sub.0.sup.h([phi])][[z([phi])[.sup.2]]/[[eta]([phi])]]dz (32)

Since the upper integration limit is a function of [phi], we employ the domain parameterization method from Ref. 28 to the fixed domain r [member of] [-1, 1] using z = [h/2](r + 1) to compute [[partial derivative]S]/[[partial derivative][[phi].sub.i]]. The design derivative of S is obtained from Eq. 32 following mathematical manipulation as

[[partial derivative]S]/[[partial derivative][[phi].sub.i]] = [[3h]/4][[[partial derivative]h]/[[partial derivative][[phi].sub.i]]]S - [[h.sup.2]/4][[integral].sub.0.sup.h][[[partial derivative][eta]]/[[partial derivative][[phi].sub.i]]][[z.sup.2]/[[eta].sup.2]]dz (33)

where the integral is evaluated with Gauss quadrature to be consistent with the numerical procedure described for Eq. 14 when evaluating S in Eq. 5. Differentiating the viscosity [eta] with respect to the design variable [[phi].sub.i] for a GNF yields

[[partial derivative][eta]]/[[partial derivative][[phi].sub.i]] = [beta][1/h][[[partial derivative]h]/[[partial derivative][[phi].sub.i]]] (34)

where the viscosity derivative factor [beta] for each GNF is given in Table 2. Note that [beta] is identical to that used in Eq. 16 when computing [[partial derivative][eta]]/[[partial derivative]P].

Therefore, [[partial derivative]R]/[[partial derivative][[phi].sub.i]] for a GNF is computed from the discretization of Eq. 29 using Eqs. 33 and 34. The design sensitivity of [g.sub.1] and [g.sub.2] in Eq. 1 then follows from Eq. 22 using Eqs. 24 and 25, respectively.

COAT HANGER DIE DESIGN EXAMPLE

To illustrate the die optimization methodology described above, we consider the design of a coat-hanger sheeting die cavity for the flow of an isothermal power-law and Carreau-Yasuda fluid. The power-law fluid is chosen here for its simplicity and since it is widely used in sheeting die design calculations. Design computations using the Carreau-Yasuda fluid model demonstrate the analysis and design sensitivity analysis given above and illustrate the effect of representing the low-strain rate near-constant viscosity in the constitutive model. Ellis and Cross fluids are not included in these example calculations since low-strain rate viscosity effects are exposed by the Carreau-Yasuda fluid analyses.

[FIGURE 2 OMITTED]

The sheeting die design considered in these examples is derived from the coat-hanger die, which is common in industrial application, and is widely studied and tested experimentally. Due to symmetry, we only mesh half of the die cavity, as illustrated in Fig. 2a. The die cavity geometry considered here follows that presented elsewhere by Gifford [10]. It consists of four major regions: the manifold, preland, secondary manifold, and the land. The purpose of the manifold is to distribute the polymer melt uniformly across the die. The preland, secondary manifold, and the land each have uniform cavity half-heights, and act as a resistance to the flow, which provides better flow uniformity [29]. The land defines the thickness of the polymer melt immediately exiting the die. The dimensions that define our die cavity geometry appear in Table 1 of Gifford [10], except for the land gap, which is fixed at 1.2 mm (i.e., [h.sub.exit] = 0.6 mm) in this study. The total die exit width is 1016 mm, which results in an exit width-to-height aspect ratio of 423. The die inlet gap and width are 19.05 and 101.6 mm, respectively, and the total die length (including the inlet channel) is 330 mm. The secondary manifold gap is fixed at 6.1 mm.

The finite element mesh with 846 nodes and 1146 elements used in this research appears in Fig. 2b. The inlet channel, secondary manifold, and land are modeled with 4-node isoparametric quadrilateral elements, while the manifold and preland contain 3-node triangular elements. Figure 2b also shows the design region, where the design variables in [phi] from Eq. 1 define the die cavity half-height h(x, y) through 445 nodal half-height design variables (bounded by 0.0001 mm [less than or equal to] [h.sup.I] [less than or equal to] 19.05 mm in the optimization). The die inlet pressure [P.sub.in] (bounded by 1.0 MPa [less than or equal to] [P.sub.in] [less than or equal to] 20 MPa in the optimization) is also defined by the design variable vector [phi], which results in a total of N = 446 design variables in the optimal design problem given in Eq. 1. The initial die cavity half-height distribution for all of the die design calculations below follows the die geometry in Gifford [10] and appears in Fig. 3a.

The objective and constraint functions are given in Eq. 1 for each design optimization considered below. Exit velocity tolerances for constraints [g.sub.1] and [g.sub.2] are defined as [[epsilon].sub.1] = 0.275 X [10.sup.-3] and [[epsilon].sub.2] = 0.225 X [10.sup.-3], respectively. In addition, [N.sub.s] = 831 height-gradient constraints are defined as ||[nabla][h.sub.p]|| = 0.25. The sequential linear programming (SLP) algorithm in Design Optimization Tools (DOT) [30] is used to solve the optimization problems.

The die cavity geometry parameterization considered in this study describes a general die gap distribution that is notably unconventional in sheeting die application. This parameterization allows for more flexibility in the design process, which is expected to yield die geometries that better satisfy Eq. 1 than if restrictions were imposed to limit die changes to those consistent with industrial applications. As a result, unconventional machining techniques may be required to fabricate the manifolds computed in the following examples. Furthermore, experiments to illustrate the accuracy of the computed die cavities given below would increase confidence in these results, but are beyond the scope of the current study.

Example 1: Die Design With an Isothermal Carreau-Yasuda Fluid

The first example considers the optimal design of the die cavity geometry described above using the Carreau-Yasuda fluid model given in Table 1. The flow of low density polyethylene (LDPE) at 270[degrees]C is selected where material constants [[eta].sub.0] = 800 Pa-s, [[eta].sub.[proportional]] = 0, [lambda] = 0.02129 s, n = 0.45958, and a = 2 are taken from Gifford [10].

Optimization results are summarized in Table 3, where the optimal die cavity design is obtained in 11 SLP design iterations. The optimization history for the inlet pressure appears in Fig. 4a and values of the constraints [g.sub.1] and [g.sub.2] in Eq. 1 are shown in Fig. 4b at each optimization iteration. In these calculations, the inlet pressure is reduced slightly by 3.8%, while the exit velocity variation constraint [g.sub.1] is reduced considerably from its initial value of 1.79 X [10.sup.-3] to its optimal value of 0.021 X [10.sup.-3], which is well below the tolerance [[epsilon].sub.1]. Similarly, the exit flow rate constraint [g.sub.2] is reduced from 2.2 X [10.sup.-3] to 0.169 X [10.sup.-3] in the same calculations. Figure 5a illustrates the gap-wise average exit velocity for the initial and optimal die designs where the distance along the die exit is normalized by dividing by the exit width of 508 mm. Note that the velocity for the optimal design is uniform over most of the die at a value that is slightly below the target velocity of 600 mm/s, with the exception of a slight reduction at the outer edge. The average exit velocity [v.sub.a], and the maximum and minimum gapwise average exit velocities are given in Table 3 to further emphasize the success of the optimization procedure.

[FIGURE 3 OMITTED]

Figure 3b gives the optimal half-height distribution in the die cavity where slight changes in the die cavity gap occur during the optimization calculations. Figure 6a shows the change in half-height h from the initial design to the optimal design and is provided here to better illustrate the changes in the design region shown in Fig. 2b. The optimal design suggests a slight increase in height along the back wall of the flow channel, which extends across the width of the die, while reducing its height downstream. An increased gap along the centerline of the die is also shown. Not as obvious from these figures is a slight decrease in the preland gap near the center of the die. The preland gap increases with x such that the optimal design has a slightly larger gap toward the outer edge of the die.

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

Example 2: Die Design With an Isothermal Power-Law Fluid

As a second example, we consider the die design problem described above using an isothermal power-law fluid. To provide a direct comparison with the Carreau-Yasuda fluid calculations above, we define m = 6402 Pa-[s.sup.n] and n = 0.45958 to match the high strain-rate viscosity in the example considered above. Results from the optimization procedure are also summarized in Table 3, where the inlet pressure is reduced by 2.9% from 10.4 to 10.1 MPa while the exit velocity constraints [g.sub.1] and [g.sub.2] are decreased below the imposed constraint tolerances given above. The change in half-height within the design region is illustrated in Fig. 6b for the power-law design. Note that changes in die cavity half-height during the power-law optimization are slightly different compared to those computed in the Carreau-Yasuda analysis above; however, the overall die cavity modifications are similar.

Exit velocity values for the power-law fluid simulations are also given in Table 3, which show improvements similar to those described above for the Carreau-Yasuda evaluations. Note that the exit velocity in the initial design is slightly higher for the power-law fluid along the die center line, and slightly higher along the die outer edge for the Carreau-Yasuda fluid, which is consistent with three-dimensional results computed by Dooley [8]. Gap-wise average exit velocities along the die exit for the power-law fluid example show similar trends as those presented for the Carreau-Yasuda fluid as illustrated in Fig. 5b.

Effect of the GNF Model

The optimization problems presented above illustrate that die cavity geometries can be defined to generate a uniform exit velocity at a desired exit velocity for various GNF fluid models. As noted above, the die cavity geometries differ only slightly between the optimal Carreau-Yasuda and optimal power-law designs. Note, however, that the optimal inlet pressure for the Carreau-Yasuda design at 9.62 MPa is 4.75% lower than that computed for the optimal power-law fluid design at 10.1 MPa. The higher inlet pressure requirement for the power-law fluid design is as expected since this fluid model overestimates the melt viscosity at low shear rates.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

To better assess the effect of the GNF model on the sheeting die design and its performance, several simulations were performed at various die flow rates. Die flow simulations were performed using the initial die cavity design shown in Fig. 3a with inlet pressure values ranging from 4 to 16 MPa in 2 MPa increments. The solid curves in Fig. 7 illustrate the computed average exit velocity [v.sub.a] for the Carreau-Yasuda and power-law fluids described above. As expected, the power-law fluid exhibits a linear pressure versus flow rate relationship on a log-log plot. Alternatively, the same plot for the Carreau-Yasuda fluid deviates from linearity with a decreasing inlet pressure as the flow rate is reduced. The difference between the power-law fluid and Carreau-Yasuda fluid response increases as the pressure and [v.sub.a] decrease. Figure 8 shows the values of [g.sub.1] from Eq. 1 for the same extrusion simulations and melt material properties. These results show that the die exit velocity uniformity remains invariant with inlet pressure (or flow rate) for the power-law fluid simulations as measured by [g.sub.1]. Alternatively, when the Carreau-Yasuda fluid is employed, the exit velocity variation decreases with inlet pressure (and flow rate). As the inlet pressure increases, values of [g.sub.1] approach those computed from the power-law simulations.

In addition to evaluating the exit velocity at different inlet pressures for a single die cavity geometry, we also consider similar computed results for die cavity designs that are optimized for each flow rate. The design optimization problem in Eq. 1 was solved several times using both the Carreau-Yasuda and power-law fluid described above for [v.sub.a] values ranging from 200 to 1000 mm/s in increments of 200 mm/s. The dashed curves in Fig. 7 illustrate the optimal inlet pressure/average velocity relationship for the die parameterization considered in this study. In this manner, we generate a set of Pareto optimal die cavity designs in a manner similar to that described by Stadler [25], which provide the least possible [P.sub.in] for a given [v.sub.a]. These Pareto optimal curves provide the design tradeoff relationship whereby the inlet pressure [P.sub.in] may only be reduced by decreasing the die flow rate (indicated here by [v.sub.a]). In other words, no die cavity geometry exists, within the design parameterization considered, that will provide a higher flow rate for the same or lower inlet pressure, while producing a uniform exit velocity. While the curve of Pareto optimal designs are useful for the die design engineer, we see from Fig. 7 that the inlet pressure is only slightly lower than that obtained from the initial design for a given [v.sub.a]. However, unlike the initial design results, the Pareto optimal designs generate a more uniform exit velocity as measured by [g.sub.1] in Eq. 1, which is illustrated by the dashed lines in Fig. 8. Die designs for both GNF fluids considered here are computed to have significantly lower [g.sub.1] values, and thus more uniform exit velocities than is seen in the initial design.

CONCLUSION

This article presents an efficient method for the design of polymer sheeting dies based on the integration of finite element flow simulations, sensitivity analysis, and optimization. Simulations that are employed in the optimal design process are based on the generalized Hele-Shaw flow model to provide a means for including die cavity geometries with arbitrary in-plane features without the computational burden characteristic of more complex three-dimensional simulations. This research extends earlier optimization-based approaches to include various GNF models, such as the power-law, Carreau-Yasuda, Cross, Ellis, and Bingham, in the analyses and sensitivity analyses required for the optimization calculations. A viscosity derivative factor is derived which facilitates both the Newton-Raphson iteration procedure and our design sensitivity computations. Optimizations are performed for power-law and Carreau-Yasuda fluids, where the die cavity geometry is based on an industrial coat hanger sheeting die design.

[FIGURE 8 OMITTED]

These computations illustrate that uniform exit velocities can be obtained via an optimization-based approach when GNF models are employed. The optimal die cavity geometries are dependent on the fluid model; however, only minor differences are seen in the final half-height values. Optimal inlet pressures are reduced when compared to the initial design, and pressures from power-law fluid optimizations exceed those computed when a Carreau-Yasuda fluid is considered, as expected. In addition, it is shown that the exit velocity uniformity depends on flow rate for a Carreau-Yasuda fluid, while it remains invariant when a power-law fluid is employed in the flow calculations. Finally, an optimization procedure is used to generate Pareto optimal designs that indicate the tradeoff between pressure drop and flow rate.

TABLE 3. Optimization results summary. Flow model, Carreau-Yasuda Power-Law Die design Initial Optimal Initial Optimal [P.sub.in] (MPa) 10.0 9.62 10.4 10.1 Constraint [g.sub.1] 1.79 0.021 3.93 0.261 (X[10.sup.-3]) Constraint [g.sub.2] 2.20 0.169 5.60 0.147 (X[10.sup.-3]) [v.sub.a] (mm/s) 628.1 592.2 644.9 592.7 [v.sub.exit] (mm/s) (max) 651.8 594.7 682.6 602.8 (min) 569.8 582.0 558.2 562.7

Contract grant sponsor: NSF; contract grant numbers: DMI-0094011; DMI-0327732.

REFERENCES

1. J.F. Carley, J. Appl. Phys., 25(9), 1118 (1954).

2. J.R.A. Pearson, Trans. J. Plast. Inst., 32, 239 (1964).

3. J.M. McKelvey and K. Ito, Polym. Eng. Sci., 11(3), 258 (1971).

4. C.I. Chung and D.T. Lohkamp, Mod. Plast., 53, 52 (1976).

5. H.H. Winter and H.G. Fritz, Polym. Eng. Sci., 26, 543 (1986).

6. J.P. Pan, P.Y. Wu, and T.J. Liu, Polym. Eng. Sci., 37(5), 856 (1997).

7. T.J. Liu, L.D. Liu, and J.D. Tsou, Polym. Eng. Sci., 34(7), 541 (1994).

8. J. Dooley, SPE ANTEC Tech. Papers, 36, 168 (1990).

9. W.A. Gifford, Polym. Eng. Sci., 40(9), 2095 (2000).

10. W.A. Gifford, Polym. Eng. Sci., 41(11), 1886 (2001).

11. Z. Tadmor and C.G. Gogos, Principles of Polymer Processing, John Wiley and Sons, New York (1979).

12. R.B. Bird, R.C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, John Wiley and Sons, New York (1987).

13. J. Nobrega, O.S. Carneiro, F.T. Pinho, and P.J. Oliveria, SPE ANTEC Tech. Papers, 48, 122 (2002).

14. J. Nobrega, F. Pinho, P. Oliveria, and O. Carneiro, SPE ANTEC Tech. Papers, 49, 310 (2003).

15. L. Sartor, "Slot Coating: Fluid Mechanics and Die Design." Ph.D. thesis, University of Minnesota (1990).

16. D.E. Smith, D.A. Tortorelli, and C.L. Tucker, SPE ANTEC Tech. Papers, 41, 42 (1995).

17. D.E. Smith, D.A. Tortorelli, and C.L. Tucker, Comput. Methods Appl. Mech. Eng., 167(3-4), 283 (1998).

18. D.E. Smith, D.A. Tortorelli, and C.L. Tucker, Comput, Methods Appl. Mech. Eng., 167(3-4), 303 (1998).

19. D.E. Smith, Int. J. Numer. Methods Eng., 57, 1381 (2003).

20. S.H. Wen and T.J. Liu, Polym. Eng. Sci., 35(9), 759 (1995).

21. C.A. Hieber and S.F. Shen, J. Non. Newtonian Fluid Mech., 7, 1 (1980).

22. P. Kennedy, Flow Analysis of Injection Molds, Hanser Publishers, New York (1995).

23. D.E. Smith, SPE ANTEC Tech. Papers, 49, 315 (2003).

24. S.J. Bates, J.F.T. Pittman, J. Sienz, and D.S. Langley, Polym. Eng. Sci., 43(8), 1500 (2003).

25. W. Stadler, J. Appl. Mech. Trans. ASME, 44, 291 (1977).

26. J.S. Arora, Introduction to Optimum Design, McGraw Hill, New York (1989).

27. K.H. Huebner, D.L. Dewhirst, D.E. Smith, and T.G. Byrom, The Finite Element Method for Engineers, John Wiley and Sons, New York (2001).

28. R.B. Haber, in Computer Aided Optimal Design: Structural and Mechanical Systems, NATO ASI Series, Vol. F27, C.A. Mota Soares, ed., Springer-Verlag, Berlin, Germany (1987).

29. D.G. Baird and D.I. Collias, Polymer Processing Principles and Design, Butterworth-Heinemann, Newton, MA (1995).

30. Vanderplaats Research and Development, Inc., Design Optimization Tools Users Manual, Vanderplaats Research and Development, Colorado Springs, CO (2004).

Douglas E. Smith, Qi Wang

University of Missouri-Columbia, E2411 EBE/MAE Department, Columbia, Missouri 65211

Correspondence to: D. Smith; e-mail: SmithDoug@Missouri.EDU

Printer friendly Cite/link Email Feedback | |

Author: | Smith, Douglas E.; Wang, Qi |
---|---|

Publication: | Polymer Engineering and Science |

Date: | Jul 1, 2005 |

Words: | 8013 |

Previous Article: | Experimental and numerical investigation of warpage of semicrystalline polymers in rotational molding. |

Next Article: | Comparison of the performance of vulcanized rubbers and elastomer/TPE/iron composites for less lethal ammunition applications. |