# A surface model-based simulation for warpage of injection-molded parts.

INTRODUCTIONInjection molding is widely used for mass producing discrete plastic parts of complex shape cost effectively with high precision. In injection molding, warpage of the molded part is possibly the most common and difficult problem. Various factors regarding the part and mold designs as well as the material selection and process setup influence the warpage of the final part. Therefore, the quality of the injection-molded product is inherently difficult to predict and/or control without employing sophisticated computer simulation software during the design stage [1].

Numerical simulation for plastic injection molding is arguably the most successful example in the field of process simulation. Considering the fact that most injection-molded components are thin walled, conventional mathematical model for injection molding simulation is based on the "mid-plane" approach [2, 3], in which the planar two-dimensional (2D) mid-plane mesh (middle-plane/center plane) with a defined thickness is used to represent the 3D geometry of the part. So, it is referred as a 2.5D model. Specific to warpage simulation, the mesh used for analysis comprises shells with reference surfaces located at the mid plane of the component.

Although the mid-plane model has been successful in predicting the fluid flow field and structural response in most cases, it does have significant implications for any useful simulation system [4]. Nowadays, part and mold designers perfect their skills at creating 3D solid geometric models using one of the popular 3D computer-aided-design (CAD) products (e.g., Unigraphics), and almost all parts are created and stored as solid-based models. So for simulation, the mid-plane model has to be extracted from the 3D solid model, or alternatively, the mid-plane mesh has to be converted from auto-meshed 3D mesh from the solid model. However, both of the tasks are not generally a simple matter and are very time consuming. In addition to the large human resources needed, potential delays in the product development process are common because this process requires continuous communication between CAD and CAE (computer-aided-engineering) models [5].

Considering the difficulty of generating a mid-plane mesh, the surface model has been introduced in injection molding simulation, in which a 3D part is represented by a boundary (skin) mesh on the outside surfaces of a solid geometric model, instead of the mid plane. Actually the surface model is still a 2.5D model, but it allows users to analyze solid geometric models of thin-walled parts directly, resulting in a significant decrease in model preparation time [6].

The surface model technique has already been used for filling, packing, and cooling simulation of the injection molding process [7-9]. For warpage, how to use the surface mesh for structural analysis is still a problem. Although Moldflow has a surface model package, it is used as a "black box" without open literatures discussing the method employed. In this article, the surface model technique is extended to the warpage simulation of injection molding, in which the molded part is decomposed into "bonded" eccentric hair-thickness shell plates. A triangular, 3-node, 6-DOF (degree of freedom) per node, flat shell element is developed for the plate structural analysis, by combining a 3-node, 3-DOF per node (including a rotational degree of freedom) membrane element and a 3-node, 3-DOF per node plate bending element. The bonding of the opposite plates is assured by multipoint constraints, and a Lagrange multiplier-based elimination method is proposed for constraint application.

BASIC IDEA

Warpage simulation of injection-molded parts is a structural analysis of shells or plates considering the accumulated residual stress. In the formulations of shell elements, there are three distinct choices: the flat shell element, the curved shell element, and the 3D solid theory-based degenerated shell element. Among them, the flat shell element is the most popular approach, as it is simplest, computationally efficient, and avoids complex shell equations.

Surface model-based warpage simulation does not need to build the mid-plane any more. For a plate shown in Fig. 1, it can be represented as a perfect bonding of two plates each with a half thickness 0.5 h. The top plate is modeled using an eccentric shell element with the top surface as its reference surface, and similarly, the bottom plate is modeled using an eccentric shell element with the bottom surface as its reference surface. To get correct structural response, these two plates are "bonded together." The bonding involves imposing the Love-Kirchh-off assumption of classical plate or shell theory, i.e., the normal to the plates remains straight after deformation and unchanged in length.

[FIGURE 1 OMITTED]

The mesh on the top surface is usually not coincident with the bottom mesh because they are generated, respectively. Hence, a normal from a mesh node Q on the bottom surface will not always find a corresponding node on the top surface. Instead, it is more likely that the normal will intersect a top element at a point P, as shown in Fig. 2. In this case, the constraints can be interpolated by using the three nodes defining the top element. To define the relationship between the degrees of freedoms of the node Q and those of its matching element, it is assumed that the normals to the reference surface before deformation remain straight after deformation. Consequently, the following relationships between the degrees of freedom of the node Q and displacements and rotations of its matching point P are adopted, as

[FIGURE 2 OMITTED]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

where z direction in the local coordinate system is the normal of the local reference surface; [u.sub.xq], [u.sub.yq], [u.sub.zq], [[theta].sub.xq], [[theta].sub.yq], and [[theta].sub.zq] are the local degrees of freedom at the node Q referred to the local element coordinate system of its matching element E; and [u.sub.xp], [u.sub.yp], [u.sub.zp], [[theta].sub.xp], [[theta].sub.yp], and [[theta].sub.zp] are the displacements and rotations at the point P in the same coordinate system.

The last relationship in the above equation defining [[theta].sub.zq] is obtained using the following drilling degree of freedom defined in the shell element formulation.

[[theta].sub.z] = [1/2][[[[partial derivative][u.sub.y]]/[[partial derivative]x]] - [[[partial derivative][u.sub.x]]/[[partial derivative]y]] (2)

Equation 1 can be rewritten in the matrix form as follows:

[U.sub.q] = [A.sub.1][U.sub.p] + [A.sub.2][D.sub.0] (3)

where

[U.sub.q] = [{[u.sub.xq], [u.sub.yq], [u.sub.zq], [[theta].sub.xq], [[theta].sub.yq], [[theta].sub.zq]}.sup.T] (4a)

[U.sub.p] = [{[u.sub.xp], [u.sub.yp], [u.sub.zp], [[theta].sub.xp], [[theta].sub.yp], [[theta].sub.zp]}.sup.T] (4b)

[D.sub.[theta]] = [{[partial derivative][[theta].sub.xp]/[partial derivative]x, [partial derivative][[theta].sub.yp]/[partial derivative]y}.sup.T] (4c)

[A.sub.1] and [A.sub.2] are matrices that depend on the coordinates of the node Q and the point P.

From the element formulation, the displacements and rotations at the point P can be obtained by shape function interpolations as follows:

[U.sub.p] = [B.sub.1][[alpha].sub.e] (5)

where [[alpha].sub.e] = [{[u.sub.xi], [u.sub.yi], [u.sub.zi], [[theta].sub.xi], [[theta].sub.yi], [[theta].sub.zi], [u.sub.xj], [u.sub.yj], [[theta].sub.xj], [[theta].sub.yj], [[theta].sub.zj], [u.sub.xk], [u.sub.yk], [u.sub.zk], [[theta].sub.xk], [[theta].sub.yk], [[theta].sub.zk]}.sup.T] and [B.sub.1] is a matrix whose values also depend on the local coordinates of the point P. Based on Eq. 2, [D.sub.[theta]] can be expressed as follows:

[D.sub.[theta]] = [B.sub.2][[alpha].sub.e] (6)

where [B.sub.2] is a matrix whose values depend on the local coordinates of the point P.

Combining Eqs. 3, 5, and 6, the relationship between the degrees of freedom of the node Q and degrees of freedoms of its matching element E is obtained as follows:

[U.sub.q] = ([A.sub.1][B.sub.1] + [A.sub.2][B.sub.2] (7)

The above constraint equations can be transformed to the global coordinate system thereby providing the final constraint equations between degrees of freedom of the node Q and degrees of freedom of its matching element E in the global coordinate system. With these constraints at all nodes on the bottom (or top) surface of the geometry, the structural performance of the composite structure can be identical to the original plate.

Based on the above modeling, the surface model-based simulation for warpage of injection-molded parts can be accomplished by use of bonding and constraints as follows: (1) building shell elements with their reference surfaces located at the outside surfaces (outer boundary) of the part; (2) using multipoint constraints to ensure that normals to the top and bottom surfaces (shell elements) remain straight after deformation.

SHELL ELEMENT FORMULATION

One approach for development of flat shell elements is to include membrane and bending properties by combining a membrane element and a plate bending element [10, 11]. Based on this approach, flat shell elements are easier to formulate using previously available theories of membrane and plate bending elements.

The constant strain triangle (CST) element [12, 13] and the linear strain triangle (LST) element [14] are early developed membrane elements. Due to the lack of drilling degrees of freedom, both of them will cause a rotational singularity in the stiffness matrix, when all the elements sharing the same node are coplanar and the local coordinate systems of the elements coincide with the global coordinate system. The first successful triangular membrane element with drilling freedoms was presented by Allman [15], followed by some similar ones [16-19]. Most of them suffer from aspect ratio locking which may occur when the response of the element is highly dependent on its aspect ratio. A direct fabrication approach was proposed by Felippa et al. to construct optimal versions of LST elements [20-22], in which the technique of Assumed Natural DEviatoric Strain (ANDES) formulation was used. Based on the ANDES template, various elements can be formulated by assigning different free parameters, e.g., an optimal membrane triangular element (OPT element) [23].

As for plate bending elements, Bazeley et al. [24] proposed a triangular plate bending element by using shape functions based on the area coordinates, called after the authors' initials as BCIZ. The nonconforming element does not pass the patch test for some mesh patterns. The Discrete Kirchoff Triangular (DKT) element is probably the most influential bending element so far [25-27], which is valid only for thin plates. The Reissner-Mindlin plate theory takes shear deformation into account in the formulation so that both thick and thin plates can be analyzed by one element model [28-31], However, there exists a shear locking phenomenon in this parametric displacement-based finite element model as the plate becomes thinner. Recently, a refined nonconforming element method (RNEM) was presented by Chen and Cheung [32, 33], based on which some refined elements have been developed and show high performance [34-39].

Based on the new progress of membrane and bending elements [23, 35], a triangular, 3-node, 6-DOF per node, i.e., three translational and three rotational degrees of freedom at each node, flat shell element is developed in this article, by combining a 3-node, 3-DOF per node (including a rotational degree of freedom) membrane element and a 3-node, 3-DOF per node plate bending element. The membrane component is derived from the ANDES template, and the employment of drilling degrees of freedom improves the accuracy of membrane response and solves the aspect ratio locking problem. By introducing the displacement function of the Timoshenko's beam into the plate bending component, the RNEM-based element possesses high accuracy for thin/thick plates.

As mentioned earlier, the flat shell element can be separated into a membrane element and a plate bending element. By assembly of the stretching-stiffness matrix and bending-stiffness matrix, the integrated stiffness matrix is as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

where [k.sub.ij.sup.m] and [k.sub.ij.sup.p] are the corresponding elements of the stretching-stiffness matrix and the bending-stiffness matrix, respectively.

The Membrane Component

Consider a flat triangular element with the corners defined by the coordinates {[x.sub.i], [y.sub.i]}, i = 1,2.3. The coordinate differences are abbreviated as [x.sub.ij] = [x.sub.i] - [x.sub.j] and [y.sub.ij] = [y.sub.i] - [y.sub.j]. The length of the side i - j is [l.sub.ij] = [square root of ([x.sub.ij.sup.2]+[y.sub.ij.sup.2])]. The thickness, area, and volume of the element are represented by h, A, and V, respectively, and [[zeta].sub.1], [[zeta].sub.2] and [[zeta].sub.3] are the area coordinates.

The degrees of freedom of the membrane component are collected in the nodal displacement vector as follows:

[{[q.sub.m]}.sup.T] = {[u.sub.1] [v.sub.1] [[theta].sub.21] [u.sub.2][u.sub.2][[theta].sub.Z2] [u.sub.3] [V.sub.3] [[theta].sub.Z3]} (9)

Based on the two-stage direct fabrication method [23], the element stiffness is the sum of the basic and the higher order stiffness matrices, as

[K.sub.m] = [K.sub.b] + [K.sub.h] (10)

where [K.sub.b] is the basic stiffness, which takes care of consistency, and [K.sub.h] is the higher order stiffness, which takes care of stability and accuracy.

An explicit form of the basic stiffness for the LST-3/9R configuration can be expressed as follows [16]:

[K.sub.b] = (1/V)[LDL.sup.T] (11)

where D is the elasticity matrix, and L is defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

where [[alpha].sub.b] is a free parameter. If [[alpha].sub.b] = 0, the basic stiffness reduces to the stiffness of the CST element.

The higher order stiffness comes from the contribution of the hierarchical corner rotations, as

[K.sub.h] = [[integral].sub.[ohm]][B.sup.T] DBd[ohm] (13)

B = [T.sub.e]([Q.sub.1][[zeta].sub.1] + [Q.sub.2][[zeta].sub.2] + [Q.sub.3][[zeta].sub.3])[~.T.sub.[theta]u] (14)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

Using Gauss integration, the final form of [K.sub.m] can be obtained as follows:

[K.sub.m] = [1/V][LDL.sup.T] + [3/4][[beta].sub.0][~.T.sub.[theta]u.sup.T][K.sub.[theta]][~.T.sub.[theta]u] (18)

where

[K.sub.[theta]] = h([Q.sub.4.sup.T][E.sub.nat][Q.sub.4] + [Q.sub.5.sup.T][E.sub.nat][Q.sub.5] + [Q.sub.6.sup.T][E.sub.nat][Q.sub.6]) (19)

[E.sub.nat] = [T.sub.e.sup.T][DT.sub.e], [Q.sub.4] = [1/2]([Q.sub.1] + [Q.sub.2]), [Q.sub.5] = [1/2] ([Q.sub.2] + [Q.sub.3]), [Q.sub.6] = [1/2]([Q.sub.3] + [Q.sub.1]) (20)

where [[beta].sub.0] is an overall scaling coefficient, and [[beta].sub.1][[beta].sub.9] are free dimensionless parameters.

An element can be generated by specifying values to the free parameters. For an isotropic material, the following values for free parameters can lead to the optimal membrane element, as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

where [upsilon] is Poisson's ratio.

The Plate Bending Component

The degrees of freedom of the plate bending component are as follows:

[{[q.sub.p]}.sup.T] = {[w.sub.1] [[theta].sub.x1] [[theta].sub.y1] [w.sub.2] [[theta].sub.x2] [[theta].sub.y2] [w.sub.3] [[theta].sub.x3] [[theta].sub.y3]} (22)

Considering the reconstruction of shear strain, its stiffness matrix includes a bending part and a shear strain part, written as follows:

[K.sub.p] = [K.sub.b] + [K.sub.s] (23)

where [K.sub.b] and [K.sub.s] are the bending part and shear part of the element stiffness matrix, respectively.

For the bending part, by eliminating the surplus parameters at the mid-side nodes as the procedure of the DKT element, the explicit expression of the rotations [[theta].sub.x] and [[theta].sub.y] for describing bending strains of the element can be obtained as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

for shape function [~.N.sub.1]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

where [[mu].sub.ij] = 1/(1 + 12[[lambda].sub.ij]), [[lambda].sub.ij] = [h.sup.2]/[5[l.sub.ij.sup.2](1-[upsilon])]; and [m.sub.ij], and [c.sub.ij] are the direction cosine and sine of the side i - j, respectively; [N.sub.i] is the shape function of the 6-node triangular element in area coordinates ([[zeta].sub.i]), as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

Other shape functions [~.N.sub.i] (i = 2, 3) can be obtained by cyclic permutation.

The bending strain of the element can be written as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

Substituting Eq. 24 into Eq. 28, the bending strain s can be expressed [epsilon] follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (29)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The part of bending of the element stiffness matrix can be written as follows:

[K.sub.b] = [[integral].sub.[ohm]] [B.sub.b.sup.T][D.sub.b][B.sub.b.sup.T][D.sub.b][B.sub.b]d[ohm] (30)

where [D.sub.b] is the flexural rigidity of the plate.

As for the shear strain part, the shear strain can be expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

To remove the shear locking, the Timoshenko's beam function is used to define the rotation and deflection on the element boundary. For instance, on the side 1-2, they read as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the tangential slope at the node i.

Thus gives as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

The displacements [~.[theta].sub.s] - [[partial derivative][bar.w.sub.s]]/[[partial derivative]l] on the side will stay constant because

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (35)

Therefore, the shear strains at node i can be expressed by the constant shear strains [gamma][S.sub.j]. For node 1, we get the following:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

Other nodal shear strains [[gamma].sub.xj] and [[gamma].sub.yj] (j = 2, 3) can be obtained by cyclic permutation similarly. Then, the shear strain can be written as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

where [bar.c.sub.i] = [c.sub.i]/[A.sub.1], [~.c.sub.i] = [c.sub.i]/[A.sub.2], [^.c.sub.i] = [c.sub.i]/[A.sub.3], [bar.m.sub.i] = [m.sub.i]/[A.sub.1], [~.m.sub.i] = [m.sub.i]/[A.sub.2], [^.m.sub.i] = [m.sub.i]/[A.sub.3], [A.sub.1] = -[c.sub.31][m.sub.12] + [c.sub.12][m.sub.31], and [A.sub.2] = -[c.sub.12][m.sub.23] + [c.sub.23][m.sub.12], [A.sub.3] = -[c.sub.23][m.sub.31] + [c.sub.31][m.sub.23], and [gamma]Sk(k = 4, 5, 6) are the natural shear strains at mid-side nodes of the element.

By further deduction, [gamma]S4 can be expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (38)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

Similarly, other mid-node parameters [[gamma].sub.S5], [[gamma].sub.S6] [B.sub.S2], and [B.sub.S3] can be obtained. The shear strains of the element can be written as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

The shear part of the stiffness matrix of the element reads as follows:

[K.sub.s] = [[integral].sub.[ohm]][B.sub.s.sup.T][D.sub.s][B.sub.s]d[ohm] (41)

where [D.sub.S] is the elasticity matrix, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], [~.D.sub.S], =5Eh/12(1+[upsilon]).

CONSTRAINT APPLICATION METHOD

In structural analysis, a constraint involving two or more nodes is called a multipoint constraint (MPC). It is obvious that the relationships expressed in Eq. 7 are typical MPCs. Generally, a MPC that connects two or more displacement components can be expressed as follows:

F([u.sub.r], [u.sub.s] ... [u.sub.t]) = prescribed value (42)

where function F involves displacement components at different nodes (two or more nodes such as r, s t, etc).

Accounting for multipoint constraints can be done by changing the assembled master stiffness equations to produce a modified system of equations, as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

The above modification process is called constraint application. The modified system is then submitted to the equation solver. The well-known methods for constraint application include the following [40]:

1. Master-Slave Elimination: The degrees of freedom involved in each MPC are separated into master and slave freedoms. The slave freedoms are then explicitly eliminated. The modified equations do not contain the slave freedoms.

2. Penalty Augmentation: Each MPC is viewed as the presence of a fictitious elastic structural clement called penalty element that enforces it approximately. This element is parametrized by a numerical weight. The exact constraint is recovered if the weight goes to infinity. The MPCs are imposed by augmenting the finite element model with the penalty elements.

3. Lagrange Multiplier Adjunction: For each MPC, an additional unknown is adjoined to the master stiffness equations. Physically, this set of unknowns represent constraint forces that would enforce the constraints exactly should they be applied to the unconstrained system.

It is emphasized that no method is uniformly satisfactory in terms of generality, robustness, numerical behavior, and simplicity of implementation. For example, the master-slave method is sensitive to the linear independence of the constraints and lacks of robustness. In cases of coupled constraint equations, incorrect decisions of slave freedoms can lead to numerical stability problems; the main disadvantage of the penalty augmentation is a serious one: the choice of weight values that balance solution accuracy with the violation of constraint conditions; the Lagrange multiplier adjunction introduces additional unknowns, requiring expansion of the original stiffness matrix and more complicated storage allocation procedures. On the whole, the third method appears to be the most elegant one for a general-purpose finite element program that is supposed to work as a "black box" by minimizing guesses and choices from its users.

The assessment summary of these MPC application methods is listed in Table 1. For a simple case with a few MPCs, these three methods are all suitable. In surface model-based warpage simulation of injection-molded parts, the part geometry is very complicated and its structural response is sensitive to the MPCs; and also, large amounts of MPCs have to be applied so as to assure the compatibility of deformation between the top and bottom plates. As a result, just using any one of the three methods cannot satisfy the requirements of warpage simulation. In this article, a Lagrange Multiplier Based Elimination (LMBE) method is presented.

TABLE 1. Assessment summary of three MPC application methods. Properties Master-slave Penalty Lagrange multiplier elimination augmentation Generality fair excellent excellent Ease of poor to fair good fair implementation Sensitivity to user high high small to none decisions Accuracy variable mediocre excellent Sensitivity to high none high constraint dependence Retains positive yes yes no definiteness Modifies unknown yes no yes vector

In the LMBE method, slave degrees of freedom are chosen for each constraint. The freedoms remaining in that constraint are labeled master. A new set of degrees of freedom ([u.sub.m]) is established only containing the constrained degrees of freedom (slave freedoms), and correspondingly, a vector [u.sub.n] is generated by removing all slave freedoms from u, i.e., it contains master freedoms as well as those that do not appear in the MPCs.

For method description, it is often convenient to transfer MPCs to the following form:

[u.sub.m] = [G.sub.mn][u.sub.n] (44)

Accordingly, the global stiffness equation can be partitioned as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (45)

To impose a multipoint constraint, imagine that the nodes involved in the constraint are connected by a rigid link, so the constraint is imposed exactly. For the convenience of computation, the link is replaced by an appropriate reaction force pair (-[lambda],+[lambda] where [lambda] is called a Lagrange multiplier. Based on the Lagrange multiplier adjunction, incorporating these forces into the original stiffness equations (as loads), and transferring it to the left hand side by appending it to the vector of unknowns, the multiplier-augmented system can be obtained as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (46)

The master stiffness matrix K is said to be bordered. Solving this system provides u and [lambda]. The latter can be interpreted as forces of constraint in the following sense: a removed constraint is replaced by a system of forces characterized by [lambda] multiplied by the constraint coefficients. More precisely, the constraint forces are -[G.sub.mn.sup.T][[lambda].sub.m].

However, it introduces additional unknowns, requiring expansion of the original stiffness matrix. Similar to the master-slave elimination, [u.sub.m] and [[lambda].sub.m] can be eliminated from the above equation, leading to the following reduced matrix as

[^.K.sub.nn][u.sub.n] = [^.f.sub.n] (47)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (48)

Equation 47 involves only the unconstrained freedoms, i.e., [u.sub.n]. The slave freedoms [u.sub.m] and additional unknowns [[lambda].sub.m] have been effectively eliminated. Solution of this modified system results in [u.sub.n], and then obtain [u.sub.m] by substituting [u.sub.n] into Eq. 44.

It should be noted that the above equation requires decisions as to which degrees of freedom are to be treated as slaves. If there are linearly dependent expressions in the constraints, [G.sub.mn] will be singular, which eventually leads to a singular [^.K.sub.nn]. To avoid this problem in the warpage simulation, we only use the nodes on the top surfaces to constrain the nodes on the bottom surfaces, that is, all constrained nodes ([u.sub.m]) are on the bottom surfaces, and all constraining nodes (in [u.sub.n]) are on the top surfaces. In addition, each node can only be constrained once. In this way, the singular problem can be avoided.

VALIDATION

First, a standard benchmark test has been carried out to evaluate the proposed model. As shown in Fig. 3, a slender cantilever beam is subjected to an end loading. The computed tip deflection (y deflection at point A) based on the current approach is compared with that of the analytical solution and ANSYS software. The results in Table 2 indicate that the proposed model maintains good accuracy for thin plates.

[FIGURE 3 OMITTED]

TABLE 2. Tip deflection for the cantilever under end loading. Method Tip deflection (mm) The proposed model 0.247 Analytical solution 0.250 ANSYS software 0.246

Then, an application in injection molding has been used to test performance of the proposed model for complicated parts. A wind board, as shown in Fig. 4a, is employed, by comparing the computation results of the proposed surface model, Moldflow's surface model and Moldflow's full 3D model. The boundary dimensions of the parts are 395.5 mm X 515.7 mm X 134.1 mm. The range of wall thickness is from 0.5 to 4.0 mm. Some ribs and bosses are designed for functional purposes. The locations signed with A, B, and C are deflection measurement positions, with the measurement datum plane determined by the locations signed with D, E, and F. The runner system and coolant layout of the mold are illustrated in Fig. 4b. The molding material used is PS-115 manufactured by Nova Chemicals. The major processing conditions are as follows: a mold temperature of 50.0[degrees]C, a melt temperature of 230.0[degrees]C, an injection time of 1.2 s, a packing profile with initial pressure of 136.0 MPa and degressive speed of 11.0 MPa/s for 10.0 s, and a cooling time of 25 s. The deflections in the experiment are measured by a vernier caliper with accuracy down to 0.01 mm.

[FIGURE 4 OMITTED]

The part deflections of numerical analyses and the experiment are illustrated in Fig. 5. It shows that the deflection trends of the three predictions are consistent with that of the experiment. The comparison of predicted and experimental deflections at the measurement points is listed in Table 3. It indicates that the deflections of the proposed surface model are smaller than the experimental data, similar to the Moldflow's surface model, whereas the deflections of Moldflow's full 3D model are larger than the experimental data. By comparison, the simulation results of the proposed model agree with the experimental data better than others.

[FIGURE 5 OMITTED]

TABLE 3. Comparison of the predicted and experimental deflections. Deflections (mm) Model Position A Position B Position C CST-BCIZ element 2.91 3.75 6.08 CST-DKT element 3.31 3.70 5.69 LST-BCIZ element 3.74 3.99 6.25 LST-DKT element 2.65 3.64 5.48 The proposed model 7.72 7.25 8.46 Moldflow's surface model 6.02 5.90 7.58 Moldflow's full 3D model 8.72 8.21 9.49 Experimental data 8.20 7.60 8.82

In addition, in order to show the performance of the presented shell element, several existing shell elements are selected for comparison with the proposed one, including the following: (1) CST-BCIZ: the 15-DOF triangular flat shell element that is derived from the membrane element CST and the plate bending element BCIZ; (2) CST-DKT: the 15-DOF triangular flat shell element that is derived from the membrane element CST and the plate bending element DKT; (3) LST-BCIZ: the 18-DOF triangular flat shell element that is derived from the membrane element LST (Allman's model with drilling degrees of freedom) and the plate bending element BCIZ; (4) LST-DKT: the 18-DOF triangular flat shell element that is derived from the membrane element LST (also Allman's model) and the plate bending element DKT. Table 3 also lists the predicted deflections of CST-BCIZ, CST-DKT, LST-BCIZ, and LST-DKT elements. The comparison shows that the presented shell element is able to produce the best results.

For comparison of computation time, the mesh size has been set as 5 mm in the proposed surface model, Moldflow's surface model, and Moldflow's full 3D model. The surface models use triangular elements, whereas the full 3D model uses tetrahedral elements. The specific element numbers of these models are shown in Table 4. The computation time for warpage analysis of the proposed surface model, Moldflow's surface model, and Moldflow's full 3D model is 1.3 min, 1.1 min, and 27 min, respectively. It can be seen that under the same element size, the element number of the full 3D model is more than 10 times that of the surface models, and accordingly, the computation time of the former is also more than 10 times that of the latter.

TABLE 4. Comparison of computation time for warpage. Element number Computation time (min) The proposed model 26164 1.3 Moldflow's surface model 25984 1.1 Moldflow's full 3D model 339096 27

CONCLUSIONS

A novel approach based on surface mesh for warpage analysis of plastic injection-molded parts has been proposed. The structure is separated into two bonded plates which are modeled as flat shell elements. The bonding of top and bottom plates is accomplished by multipoint constraints and a Lagrange multiplier based elimination method has been presented for dealing with the constraints.

The plastic injection-molded part often exhibits complex geometry (many ribs, gussets, and bosses are often designed for specific purpose). By comparing the simulation results to those of the experiment and Moldflow, it can be found that the proposed model is capable of performing a warpage analysis using only a surface mesh for injection-molded parts.

The formulation of the presented shell element is relatively simple because it is developed by superimposing the stiffness of membrane and bending elements. Another reason is that only physical degrees of freedom and corner nodes are used. The development meets the back trend toward simplicity. By comparison with some popular shell elements, the presented element shows a superior property.

REFERENCES

(1.) D.S. Choi and Y.T. Im, Compos. Struct., 47, 655 (1999).

(2.) V.W. Wang, C.A. Hieber, and K.K. Wang, J. Polym. Eng., 7, 21 (1986).

(3.) H.H. Chiang, C.A. Hieber, and K.K. Wang, Polym. Eng. Sci., 31, 116(1991).

(4.) Z. Fan, R. Zheng. P. Kennedy. H. Yu, A. Bakharev, and S.P.E. Spe, "Warpage Analysis of Solid Geometry," in Antec 2000: Society of Plastics Engineers Technical Papers, Conference Proceedings, Orlando (2000).

(5.) W.Y. Cho, J.K. Nam, S.M. Yun, and H.C. Sin, Proc. Inst. Mech. Eng. Part B: J. Eng. Manuf., 221, 845 (2007).

(6.) D. Cardozo, J. Reinf. Plast. Compos., 27, 1963 (2008).

(7.) W. Cao, R. Wang, and C.Y. Shen, Polym. Plast. Technol. Eng., 43, 1471 (2004).

(8.) H.M. Zhou and D.Q. Li, Adv. Polym. Tech., 20, 125 (2001).

(9.) H. Yu, C. Kietzmann, P. Cook, S. Xu, F. Costa, and P. Kennedy, "A New Method for Simulation of Injection Molding." in Antec 2004: Society of Plastics Engineers Technical Papers, Conference Proceedings, Chicago (2004).

(10.) E. Gal and R. Levy, Arch. Comput. Method Eng., 13, 331 (2006).

(11.) Y.X. Zhang and K.S. Kim, Comput. Mech., 36, 331 (2005).

(12.) R.W. Clough, The Finite Element Method--A Personal View of Its Original Formulation, in From Finite Elements to the Troll Platform--The Ivar Holand 70th Anniversary Volume, Tapir Academic Press, Trondheim (1994).

(13.) M.J. Turner, R.W. Clough, H.C. Martin, and L.J. Topp, J. Aerosol. Sci., 23, 805 (1956).

(14.) O.C. Zienkiewicz, Int. J. Numer. Methods Eng., 52, 287 (2001).

(15.) D.J. Allman. Comput. Struct., 19, 1 (1984).

(16.) P.G. Bergan and C.A. Felippa. Comput. Methods Appl. Mech. Eng., 50, 25 (1985).

(17.) D.J. Allman, Int. J. Numer. Methods Eng., 26, 2645 (1988).

(18.) C.C. Rankin, F.A. Brogan, W.A. Loden, and H. Cabiness, STAGS User Manual, Lockheed Mechanics, Materials and Structures Report P032594 (1998).

(19.) H.T.Y. Yang. S. Saigal, A. Masud, and R.K. Kapania, Int. J. Numer. Methods Eng., 47, 101 (2000).

(20.) K. Alvin, H.M. de la Fuente, B. Haugen, and C.A. Felippa, Finite Elem. Anal. Des., 12, 163 (1992).

(21.) C.A. Felippa and C. Militello, Finite Elem. Anal. Des., 12, 189 (1992).

(22.) C.A. Felippa and S. Alexander, Finite Elem. Anal. Des., 12, 203 (1992).

(23.) C.A. Felippa, Comput. Methods Appl. Mech. Eng., 192, 2125 (2003).

(24.) G.P. Bazeley, Y.K. Cheung, B.M. Irons, and O.C. Zienkiewiez, "Triangular Elements in Plate Bending, Confirming and Non-Confirming Solutions," in Proceedings 1st Conference on Matrix Methods in Structural Mechanics. Dayton (1966).

(25.) J.L. Meek and H.S. Tan, Comput. Struct., 21, 1197 (1985).

(26.) G. Dhatt, L. Marcotte, and Y. Matte, Int. J. Numer. Methods Eng., 23, 453 (1986).

(27.) I. Katili, Int.J. Numer. Methods Eng., 36, 1859 (1993).

(28.) J.L. Batoz and P. Lardeur, Int. J. Numer. Methods Eng., 28, 533 (1989).

(29.) J.L. Batoz and I. Katili, Int. J. Numer. Methods Eng., 35, 1603 (1992).

(30.) R.H. MacNeal, Finite Elem. Anal. Des., 30, 175 (1998).

(31.) H. Sofuoglu and H. Gedikli, Commun. Numer. Methods Eng., 23. 385 (2007).

(32.) W.J. Chen and Y.K. Cheung, Int. J. Numer. Methods Eng., 40, 3919 (1997).

(33.) W.J. Chen and Y.K. Cheung, Int. J. Numer. Methods Eng., 41, 1507 1998).

(34.) W.J. Chen and Y.K. Cheung, Int. J. Numer. Methods Eng., 51, 1259 (2001).

(35.) G. Zengjie and C. Wanji, Comput. Mech., 33, 52 (2003).

(36.) W.J. Chen, Int. J. Numer. Methods Eng., 60, 1817 (2004).

(37.) Z. Wu, R.G. Chen, and W.J. Chen, Compos. Struct., 70, 135 (2005).

(38.) Z. Wu, and W.J. Chen, Int. J. Numer. Methods Eng.. 69, 1627 (2007).

(39.) Z. Wu, W.J. Chen, and X.H. Ren, Compos. Struct., 88, 643 (2009).

(40.) O.C. Zienkiewicz, R.L. Taylor, and J.Z. Zhu, The Finite Element Method: Solid Mechanics, Butterworth-Heinemann, Oxford (2005).

Huamin Zhou, Zhiyong Wang, Jianhui Li, Dequn Li

State Key Laboratory of Materials Processing and Mold Technology Huazhong University of Science and Technology, Wuhan 430074, People's Republic of China

Correspondence to: H. Zhou; e-mail: hmzhou@hust.edu.cn

Contract grant sponsor: National Natural Science Foundation Council of China; contract grant number: 50875095; contract grant sponsor: National High-tech R&D Program (863 Program) of China; contract grant number: 2009AA03Z104.

DOI 10.1002/pen.21884

Published online in Wiley Online Library (wileyonlinelibrary.com).

[C] 2010 Society of Plastics Engineers

Printer friendly Cite/link Email Feedback | |

Author: | Zhou, Huamin; Wang, Zhiyong; Li, Jianhui; Li, Dequn |
---|---|

Publication: | Polymer Engineering and Science |

Article Type: | Report |

Geographic Code: | 9CHIN |

Date: | Apr 1, 2011 |

Words: | 6189 |

Previous Article: | Synthesis of Y-shaped Poly(N,N-dimethylamino-2-ethyl methacrylate) and Poly(trimethylene carbonate) from a new heterofunctional initiator. |

Next Article: | Thermomechanical and surface properties of novel poly(ether urethane)/polyhedral oligomeric silsesquioxane nanohybrid elastomers. |

Topics: |