The multi-point constraint ensuring convergence of surface models for warpage analysis of injection-molded plastics

Abstract

The surface models have been widely used to predict warpage of plastics for several years, but very less attention has been paid to the convergence of this kind of models. In this article, the multi-point constraint ensuring convergence of surface models is deduced based on the theoretical basis of domain decomposition method. And a new corresponding surface model is developed to verify its effectiveness in practice. Three numerical tests are used to evaluate the performances of the new surface model. In the former two cases, the proposed surface model can rapidly converge to the exact solutions, as same as the mid-plane model. And in the last case, the results from the proposed models are much closer to the experiment data than the existing models. It proves that the proposed multi-point constraint is sufficient to ensuring convergence of surface models.

Introduction

Injection molding is a common manufacturing process to produce large numbers of plastic parts of complex shapes. ¹ Warpage is a major problem both at the design stage and during the mass production, ² and simulation tools for injection molding have been widely used to help engineers to handle it. ³ The models used in warpage prediction are the key to obtain high-precision simulation results.

For warpage analysis of injection-molded plastics, the general models can be divided into three categories: the mid-plane models, the surface models (also called dual-domain models) and the solid models. ⁴ The earliest used models for injection molding are the mid-plane models,^5,6 which are suitable for most thin-walled plastics. For complex structures, the difficulty of geometry remodeling inhibits its wide application. To avoid remodeling, the surface models have been formulated, in which a mesh on the outer boundary surfaces is used to represent three-dimensional (3D) parts. At the present time, the 3D solid models are developing fast.^7,8 However, their computation time is much longer than the surface models, ⁴ while the accuracy for thin-walled parts is not obviously improved when compared with the surface models. Therefore, the surface model is still an effective and efficient method for large and complicated thin-walled parts.

In surface models, a plate is divided into two identical subdomains along the thickness, first, each represented by a shell. Then, the two subdomains are bonded by enforcing multi-point constraint (MPC). From the geometric point of view, the plate can be seen as perfect bonding of two matching shells; from the mathematical point of view, the above-mentioned process can be considered as a domain decomposition method.

There are two representative surface models for warpage analysis, proposed by Fan et al.^9,10 and Zhou and colleagues,^11,12 respectively. Both the above-mentioned models have the same MPC equation, which is based on the Love–Kirchhoff assumption. It requires that a normal to the part remains straight after deformation and the length is unchanged. But the applied shells of the two models take account of shear factor, leading to a change in the length of the normal. Hence, contradiction occurs and the convergences of the existing models would not like to be ensured certainly in all cases.

For the domain decomposition method, the convergence condition has been studied, which also has been successfully applied in the finite element method.^13,14 In the existing domain decomposition models, the interfaces between subdomains are usually of overlapping but non-matching meshes. But for a surface model, there are no meshes at the interface of two subdomains. So, the convergence condition for this model should be re-deduced from the basic MPC condition.

In this article, the MPC ensuring convergence of the surface models is derived based on the theory of the domain decomposition method. In order to validate its effectiveness, a new surface model is formulated and two benchmark cases and a practical example are implemented. The results demonstrate that the new model using the proposed MPC can always show superior behavior.

The convergence condition

The warpage problem for injection-molded parts can be described by the elliptic partial differential equation

$$\begin{array}{l}D\Delta u=f\mathrm{in}\Omega \\ u=0\mathrm{on}\partial \Omega \end{array}$$

(1)

where $u$ is the displacement, $D$ is the elastic modulus, f is force per unit volume and $\Omega$ is the domain of parts with $\partial \Omega$ as its boundary. By generating surface meshes on the outer boundary of the part, $\Omega$ is divided into two non-overlapping subdomains ${\Omega }_1$ and ${\Omega }_2$ by the mid-plane $\Gamma$ , as shown in Figure 1. Because most plastics are thin-walled, the two subdomains can be represented by two eccentric shells independently, ignoring boundary fraction. Hence, the problem described in equation (1) becomes

$$\begin{array}{cc}\multicolumn{1}{c}{}& D\Delta u_1=f\mathrm{in}{\Omega }_1\hfill \\ \multicolumn{1}{c}{}& u_1=0\mathrm{on}\partial {\Omega }_1\hfill \\ \multicolumn{1}{c}{}& D\Delta u_2=f\mathrm{in}{\Omega }_2\hfill \\ \multicolumn{1}{c}{}& u_2=0\mathrm{on}\partial {\Omega }_2\hfill \end{array}$$

(2)

where $u_i$ is the restriction of u to ${\Omega }_i(i=1,2)$ . Based on theory of the domain decomposition method, the solution of equation (2) is equivalent to that of equation (1) under the following condition ¹³

$$\begin{array}{cc}\multicolumn{1}{c}{}& u_1=u_2\mathrm{on}\Gamma \hfill \\ \multicolumn{1}{c}{}& \frac{\partial u_1}{\partial n_1}=-\frac{\partial u_2}{\partial n_2}\mathrm{on}\Gamma \hfill \end{array}$$

(3)

where $n_i$ is the outward normal to ${\Omega }_i$ in $\Gamma$ . Hence, the surface model would give equivalent solution to the corresponding mid-plane model, if equation (3) can be satisfied explicitly.

Figure 1.

A plate of a part is divided into two subdomains, each represented by a pair of eccentric shells.

From the constraint variational principle, the energy functional corresponding to the surface models can be obtained using the Ritz method and Lagrange multiplier method, ¹⁵ giving

$$\begin{array}{cc}\multicolumn{1}{c}{}& \Pi (u,{\lambda }_1,{\lambda }_2)=\frac{1}{2}\underset{\Omega }{\int }\nabla u^{T}D\nabla \mathrm{udv}-\underset{\Omega }{\int }\mathrm{fudv}\hfill \\ \multicolumn{1}{c}{}& +\underset{\Gamma }{\int }{\lambda }_1(u{|}_{{\Omega }_1}-u{|}_{{\Omega }_2})d\Gamma +\underset{\Gamma }{\int }{\lambda }_2(\frac{\partial u}{\partial n_1}{|}_{\Omega 1}+\frac{\partial u}{\partial n_2}{|}_{\Omega 2})d\Gamma \hfill \end{array}$$

(4)

With the principle of minimum potential energy, the value of the first-order variation of the functional $\Pi (u,{\lambda }_1,{\lambda }_2)$ should be equal to 0, ¹⁵ and the following two equations can be derived as

$$\underset{\Omega }{\int }\nabla \delta u^{T}D\nabla \mathrm{udv}-\underset{\Omega }{\int }f\delta \mathrm{udv}=0$$

(5)

$$\begin{array}{cc}\multicolumn{1}{c}{}& \underset{\Gamma }{\int }\delta {\lambda }_1(u{|}_{{\Omega }_1}-u{|}_{{\Omega }_2})d\Gamma =0\hfill \\ \multicolumn{1}{c}{}& \underset{\Gamma }{\int }\delta {\lambda }_2(\frac{\partial u}{\partial n_1}{|}_{\Omega 1}+\frac{\partial u}{\partial n_2}{|}_{\Omega 2})d\Gamma =0\hfill \end{array}$$

(6)

For a surface model, its stiffness matrix and load vector can be built for both the subdomains from equation (5), while equation (6) is the condition in the weak form ensuring its equivalence to the original problem.

The development of a surface model contains two parts: the flat eccentric shell representing the two separate subdomains and the MPC bonding the two domains. In this article, the triangular flat eccentric shell is adopted. It can be divided into eccentric plate component and eccentric membrane component.

Taking the “Top shell” shown in Figure 1 as an example, the reference plane (z = 0) of an eccentric shell lies in the outer boundary of the part, opposite to its neutral plane (z = −0.5 h). The displacements of the eccentric shell at point A (x, y, z) in the local coordinate are described as

$$\begin{array}{cc}\multicolumn{1}{c}{}& u{(x,y,z)}_x=u_{x0}(x,y)+(z+0.5h){\theta }_{y0}(x,y)\hfill \\ \multicolumn{1}{c}{}& u{(x,y,z)}_y=u_{y0}(x,y)-(z+0.5h){\theta }_{x0}(x,y)\hfill \\ \multicolumn{1}{c}{}& u{(x,y,z)}_z=u_{z0}(x,y)\hfill \\ \multicolumn{1}{c}{}& \theta {(x,y,z)}_x={\theta }_{x0}(x,y)\hfill \\ \multicolumn{1}{c}{}& \theta {(x,y,z)}_y={\theta }_{y0}(x,y)\hfill \\ \multicolumn{1}{c}{}& \theta {(x,y,z)}_z={\theta }_{z0}(x,y)\hfill \end{array}$$

(7)

where $u_{x0},u_{y0},u_{z0}{\theta }_{\mathrm{xo}},{\theta }_{\mathrm{yo}},{\theta }_{\mathrm{zo}}$ are the displacements of the A’s matching point in the reference plane (z = 0). And $u_{x0},u_{y0},u_{z0}{\theta }_{\mathrm{xo}},{\theta }_{\mathrm{yo}},{\theta }_{\mathrm{zo}}$ can be expressed by nodal freedoms and shape functions as

$$\left[\begin{array}{c}\hfill u_{x0}\hfill \\ \multicolumn{1}{c}{u_{y0}}\\ \multicolumn{1}{c}{u_{z0}}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{xo}}}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{yo}}}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{zo}}}\end{array}\right]=T_{6*18}{u}_{\mathrm{nodal}\mathrm{freedoms}}$$

(8)

At the element level, for two matching eccentric shells in the global coordinate as shown in Figure 2, substituting equation (7) into equation (6) gives

$$\begin{array}{c}\multicolumn{1}{c}{\int \underset{e}{\int }(u_{x0}(x,y){|}_{\mathrm{Top}}-u_{x0}(x,y){|}_{B\mathrm{ottom}})\mathrm{dxdy}=0}\\ \multicolumn{1}{c}{\int \underset{e}{\int }(u_{y0}(x,y){|}_{\mathrm{Top}}-u_{y0}(x,y){|}_{\mathrm{Bottom}})\mathrm{dxdy}=0}\\ \multicolumn{1}{c}{\int \underset{e}{\int }(u_{z0}(x,y){|}_{\mathrm{Top}}-u_{z0}(x,y){|}_{\mathrm{Bottom}})\mathrm{dxdy}=0}\\ \multicolumn{1}{c}{\int \underset{e}{\int }({\theta }_{x0}(x,y){|}_{\mathrm{Top}}-{\theta }_{x0}(x,y){|}_{\mathrm{Bottom}})\mathrm{dxdy}=0}\\ \multicolumn{1}{c}{\int \underset{e}{\int }({\theta }_{y0}(x,y){|}_{\mathrm{Top}}-{\theta }_{y0}(x,y){|}_{\mathrm{Bottom}})\mathrm{dxdy}=0}\\ \multicolumn{1}{c}{\int \underset{e}{\int }({\theta }_{y0}(x,y){|}_{\mathrm{Top}}-{\theta }_{z0}(x,y){|}_{\mathrm{Bottom}})\mathrm{dxdy}=0}& \hfill & \hfill \end{array}$$

(9)

where subscripts “Top” and “Bottom” denote the corresponding shells. With equation (8) and (9) the following relationship between two shells can be obtained as

$$\int \underset{e}{\int }{T}_{6*18}(u_{\mathrm{nodal}\mathrm{freedoms}}{|}_{\mathrm{Top}}-u_{\mathrm{nodal}\mathrm{freedoms}}{|}_{B\mathrm{ottom}})\mathrm{dxdy}=0$$

(10)

Figure 2.

The “matching nodes” N and M in a pair of opposite shells, respectively.

Based on equation (10), the following sufficient condition is proposed

$$u_{\mathrm{nodal}\mathrm{freedoms}}{|}_{\mathrm{Top}}-u_{\mathrm{nodal}\mathrm{freedoms}}{|}_{B\mathrm{ottom}}=0$$

(11)

Because the constraint in equation (11) is stronger than that in equation (6), equation (11) is the MPC ensuring convergence written in the form of nodal freedoms.

A new surface model

Based on the deduced MPC above, a new surface model is proposed. The shape functions of the outstanding membrane element OPT (optimal membrane element) ¹⁶ and the famous plate element DKT (discrete Kirchhoff element) ¹⁷ are adopted to develop the corresponding eccentric shell. With above-developed eccentric shell, the finite element equilibrium equation for two subdomains based on the equation (5) is built

$$\left[\begin{array}{cc}\hfill K_{11}\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill K_{22}\hfill \end{array}\right]\left[\begin{array}{c}\hfill u_1\hfill \\ \multicolumn{1}{c}{u_2}\end{array}\right]=\left[\begin{array}{c}\hfill f_1\hfill \\ \multicolumn{1}{c}{f_2}\end{array}\right]$$

(12)

And for a pair of matching nodes shown in Figure 2, the proposed constraint condition based on equation (11) can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{u_{\mathrm{Mx}0}(x,y)=u_{\mathrm{Nx}0}(x,y)}\\ \multicolumn{1}{c}{u_{\mathrm{My}0}(x,y)=u_{\mathrm{Nx}0}(x,y)}\\ \multicolumn{1}{c}{u_{\mathrm{Mz}0}(x,y)=u_{\mathrm{Nx}0}(x,y)}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{Mx}0}(x,y)={\theta }_{\mathrm{Nx}0}(x,y)}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{My}0}(x,y)={\theta }_{\mathrm{Nx}0}(x,y)}\\ \multicolumn{1}{c}{{\theta }_{\mathrm{Mz}0}(x,y)={\theta }_{\mathrm{Nx}0}(x,y)}\end{array}$$

(13)

Hence, for matching meshes, the constraint expressed in equation (13) can ensure convergence of the surface models certainly. And for non-matching meshes, as mesh quantity increases, the deviation would approximate to zero and the convergence still could be obtained.

The final MPC equation defined by equation (13) between subdomains ${\Omega }_1$ and ${\Omega }_2$ can be expressed as

$$A_1{u}_1=A_2{u}_2$$

(14)

where $A_1$ and $A_2$ are coefficient matrices here. Enforcing equation (14) into equation (12) with Lagrange multiplier method, the final solvable equation is

$$\left[\begin{array}{ccc}\hfill K_{11}\hfill & \hfill 0\hfill & \hfill A_1^T\hfill \\ \multicolumn{1}{c}{0}& \hfill K_{22}\hfill & \hfill -A_2^T\hfill \\ \multicolumn{1}{c}{A_1}& \hfill -A_2\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill u_1\hfill \\ \multicolumn{1}{c}{u_2}\\ \multicolumn{1}{c}{\lambda }\end{array}\right]=\left[\begin{array}{c}\hfill f_1\hfill \\ \multicolumn{1}{c}{f_2}\\ \multicolumn{1}{c}{0}\end{array}\right]$$

(15)

Eliminating three translational and three rotational degrees of freedom from the above equation, the unique solution can be obtained.

Verification

In order to evaluate the presented MPC and the corresponding surface model, two benchmark tests have been implemented to verify its responses under the membrane-dominated and bending-dominated conditions, respectively. For comparison, a mid-plane model using the shell of the same functions as the new surface model in the reference plane has been developed. In addition, the new model has also been used to predict the warpage of a practical injection-molded “wind board,” by comparing with Moldflow and experiment.

Rhombic cantilever

A rhombic cantilever subjected to a uniformed load is shown in Figure 3, whose geometry and properties are also presented. The result from Batoz et al. ¹⁷ is used for comparison. A 4 × 4 matching mesh is used for two models. The results are given in Table 1. The surface model gives nearly the same results as the mid-plane model, which proves that the proposed MPC can ensure convergence in membrane-dominated condition.

Figure 3.

Rhombic cantilever: q = 0.26066 psi, E = 10.5 × 10⁶ psi, h = 0.125 in and v = 0.3.

Table 1.

The results of two models for rhombic cantilever.

Models	Points
	1	2	3	4	5	6
The proposed surface model	0.3037	0.1986	0.1128	0.1212	0.0556	0.0225
The mid-plane model	0.3037	0.1986	0.1128	0.1212	0.0556	0.0225
Exact	0.297	0.204	0.121	0.129	0.056	0.022

Scordelis–Lo’s cylindrical roof test

The Scordelis–Lo’s cylindrical roof supported by the rigid diaphragms at two ends is shown in Figure 4. And it is severed here to assess the surface model’s ability of approximating complicated states of membrane strain. ¹⁸ The roof is subjected to a uniform gravity load of 90.0 per unit area. And the material modulus is E = 4.32 × 10⁸ and Poisson ratio is v = 0.0. The results in Table 2 show that the surface model could obtain nearly the same solutions with the mid-plane model. As mesh number increases, the results of both models tend to the exact solutions. It verifies that the proposed MPC can ensure convergence in bending-dominated condition.

Figure 4.

Scordelis–Lo’s cylindrical roof under uniform load.

Table 2.

The results of two models for Scordelis–Lo’s cylindrical roof.

Models	Mesh
	2*2	4*4	8*8	16*16
The proposed surface model	0.4544	0.3245	0.3032	0.2980
The mid-plane model	0.6104	0.3271	0.3040	0.2981
Exact	0.3024

An injection-molded wind board

An injection-molded wind board is chosen to evaluate the performance of the surface model for complicated parts, as depicted in Figure 5(a). Its geometry dimension is 395.5 mm × 515.7 mm × 134.1 mm with thickness ranging from 0.5 to 4.0 mm. The points marked with “A” and “B” are deflection-measured positions of the part, while the other three points marked as “C,”“D” and “E” are used to determine the datum plane for measurement. The major processing parameters are displayed in Table 3. The initial residual stress for load in the analysis is obtained by software system labeled as HsCAE7.5 developed by Huazhong University of Science and Technology.

Figure 5.

The “wind board”: (a) geometry, runner system and cooling system; (b) the predicted warpage field and (c) experimental deflection field.

Table 3.

The major processing parameters for molding wind board.

Material	PS-115 produced by Nova Chemicals
Mold temperature	50 °C
Melt temperature	230 °C
Injection time	1.2 s
Packing time	10.0 s
Packing pressure	Initial value: 107 MPa; degressive speed: 10 MPa/s

Figure 5(b) and (c) depicts the deflection field of the numerical prediction and experiment, and the tendencies are the same. The two points’ predicted values are compared with the measured ones and those of Moldflow, as listed in Table 4. Obviously, the surface model gives much closer results to the experimental ones than those of Moldflow. The maximum deviation percentage for all points is less than 4%, which is an acceptable result for injection warpage analysis.

Table 4.

Comparison of the experimental and computed data.

	The maximum deflection	Point defections (mm)
		A	C
The proposed surface model	8.59	7.93	8.59
Moldflow surface model	7.58	6.02	7.58
Moldflow 3D model	9.49	8.72	9.49
Experimental results	8.82	8.20	8.82

3D: three-dimensional.

Conclusion

In the article, the MPC ensuring convergence for surface models is deduced from domain decomposition method on the basis of mathematical theory. And a surface model using the proposed MPC is developed to verify its effectiveness in practice. And two benchmark tests and an injection-molded wind board are carried out to validate the performance of the proposed surface model. Compared with the mid-plane model, the surface model is of extremely small error, nearly of the same results in the benchmark tests. For the wind board example, the predicted displacements of the new model agree with the experimental data well, showing superiority to other models. Hence, the MPC is effective and the surface model using it would give convergent solutions. It is valuable for optimizing the manufacturing process of injection molding. Although the MPC in this article is proposed to ensure accurate prediction of warpage of injection-molded parts, it is not just limited to it. In other words, it can also be used to analyze general elastic deformation problems of thin-walled structures. Further work would be to determine the error variation as mesh number increases for the surface model with the MPC.