Thermoelasticity of Injection-Molded Parts

Abstract

In the realm of injection-molded parts, small length scale deformation defects such as sink marks often pose a major challenge to the aesthetics or functionality of the parts. To address this problem, we present a comprehensive thermoelastomechanical approach that calculates the deformation of injection molded plastic by solving the elastic problem at each time step. In our methodology, two treatments of the molten core are considered: one as a liquid and the other as a rubbery state. Our results suggest that the rubbery state treatment provides higher accuracy in predicting the deformation results, as it maintains the displacement of the localized thermal shrinkage in its vicinity. The validity of our method is supported by empirical measurements on produced parts from the existing literature as well as on samples that we molded independently.

1. Introduction

Injection molding of thermoplastic polymers is a widely used production process in which molten plastic is injected at high temperature into a metal cavity and cooled until it solidifies into a functional part. This process is essential to the mass production industry due to its scalability and low unit price, which keeps production costs manageable. Considering the low profit margins and high capital intensity associated with acquiring steel molds for plastic injection molding, it is critical to avoid costly mold and part design mistakes. Broken molds or visually deficient parts can have a devastating impact on profitability, which has led to the emergence of slogans such as “right on the first try”. To consistently achieve such production goals, simulation of the injection molding process is required. This involves solving the relevant partial differential equations for a given part geometry, considering appropriate boundary conditions as well as material models, and generating results that predict part moldability and quality. The prediction of warpage, i.e., the deviation of the shape of the produced part from the intended shape, is the pinnacle of modern simulation software [1,2,3].

In addition to warpage, there are many other defects that can occur in injection molded parts, and the goal of simulation is to predict as many defects as possible. One particular defect that occurs in injection molding is a “sink mark”, an undesirable depression on the surface of the part caused by localized high shrinkage of the material. Typically, sink marks occur in areas of the part that are thicker and cool more slowly, causing the surface to shrink away from the cavity wall more than adjacent areas. The local thickness is best defined as the diameter of a maximum sphere contained locally in the part, as shown in Figure 1 [4].

The aim of this work is to show that both global warpage and local depressions are due to the same cause, namely, thermal contraction during cooling and solidification of the plastic under pressure. Although this interpretation has been suspected in the literature, it has not yet been clearly established. The exact elastomechanics of the solidified layer in at least two dimensions throughout the production cycle must be calculated to show that sink marks and warpage are the same deformation manifested at different length scales. Furthermore, this work aims to broaden the definition of “sink marks” to include both depressions under a typical rib, as shown on the far left in Figure 1, and local surface deformations caused by any local variation of thickness.

1.1. Historical Overview of Injection Molding Process Modeling

An excellent historical overview of the academic field of injection molding simulation is provided by Peters Kennedy’s book [5], which serves as a common entry point for this field of study. The beginnings date back to 1960 when electronic computers were introduced and the first non-isothermal calculation of polymer flow in a simple rectangular plate cavity was performed [6]. This calculation involved modeling the increase in viscosity due to heat loss and discussing the resulting influence on the molding process. It then took 18 years of academic demonstrations before the first commercial company, Moldflow, was founded by Colin Austin in 1978. Engineers now had access to academic-grade code and the book Moldflow Design Principles [7], which still serves as the benchmark for the field. At that time, numerical techniques were used to reduce three-dimensional sheet-like geometries to equivalent flow in a plane [8].

With the increasing popularity of computer-aided design systems (CAD) in the 1980s, midplane shell models emerged. It became possible to solve the relevant equations using the finite difference method over the thickness of the part and to propagate them over the entire component using the finite element method [9]. It should be noted, however, that this type of modeling cannot capture the local variation in thickness. Academic work paved the way for modeling the packing phase, i.e., the time after the cavity is completely filled but the plastic is still molten and under pressure [10,11,12,13]. During this period, commercial companies began to develop the calculation of heat transfer both through the mold and through warpage. Moldflow led the way in warpage prediction, with a phenomenological approach, molding and simulating 28 sample plates with different parameters. The shrinkage of the sample plates was measured and used to calibrate the parameters of the “Residual Strain” model [14].

In the 1990s, CAD systems evolved to produce three-dimensional models. Full three-dimensional simulation using the “volume of fluid” method gained traction [15], but computer power remained limited and workarounds were developed. The most renowned approach is arguably dual domain technology [16], where the three-dimensional CAD geometry surface is meshed and computational cells are formed from matching mesh elements on the opposite surfaces of the part. As computing power increased, full three-dimensional simulation became viable [17], leading to a major expansion of simulation functionality. Increased accuracy enabled direct calculation of phenomena such as warpage [18] and three-dimensional velocity fields.

Open source software for injection molding simulation was virtually nonexistent until the early stages of development were begun in [19]. The resulting hydrodynamic flow simulation code [19] is based on a solver within the popular open source finite volume library OpenFOAM [20]. Over time, functionalities for the calculation of fiber orientation and crystallization were added [21]. Validation was performed [22] and good performance was shown.

1.2. Modeling the Deformation of Injection Molded Parts

During the injection molding process, the plastic is exposed to significant temperature and pressure changes as the phase boundary between solid and molten plastic progresses. Already-solidified layers prevent newly formed layers from freely changing shape, as they are all subject to elastomechanics as an integral domain. Early attempts to model shrinkage of injection molded parts relied only on thermodynamic (pvT) characterization of density [23] to predict both isotropic stress and shrinkage. Later, a line of research was devoted to viscoelastic formulations of the problem, which is complex from both theoretical and experimental point of view, as it is difficult to measure viscoelastic material properties in a satisfactory way [24]. Nevertheless, Douven [25] has shown that the stresses calculated by elastic and viscoelastic descriptions are similar. Moreover, it was shown in parallel by Baaijens [26] that the flow-induced stresses are two orders of magnitude smaller than the thermal and pressure-induced stresses. Therefore, the elastic thermomechanical approach proved sufficient, which provided the basis for its validity for a group of authors who worked on its increasingly rigorous formulation [27,28]. These efforts culminated in a paper by Jansen [29] in which a fully elastic thermomechanical method for calculating residual stress and shrinkage was developed. Their method showed that the boundary conditions at the interface between the plastic and the mold play a more decisive role than the increasingly complex material characterization. The objective of the present study is to take up Jansen’s approach and adapt it for a two-dimensional cross-section of an injection molded part. Specifically, the aim is to predict the deformation of the part surface with a characteristic length scale of the part thickness.

1.3. Modeling the Formation of Sink Marks

Sink marks were first discussed by Marchewka in 1974 [30]. In their study, plates with ribs were molded under different processing conditions and sink marks were visually assessed. The study focused primarily on the visibility of sink marks as a function of rib thickness, distance from the gate, and surface finish of the affected area. The results of the study led to the establishment of a design guideline that a rib thickness of less than $60\%$ of the base plate should not produce visible sink marks under most processing conditions. In addition, the study found that the distance from the gate is an important factor affecting the visibility of sink marks.

Later, Naka et al. [31] of what is now Phillips Corporation developed a method that used finite differences to calculate the evolving temperature distribution in a two-dimensional part cross section during injection molding. The temperature change was used as a thermal load for one-dimensional fiber elements that were independent of each other and constrained in the longitudinal direction. Only after the pressure in the cavity reached atmospheric values were the thermal, plastic, and elastic strains calculated over time and finally converted into displacements in the cross-sectional plane using the Poisson effect for the elastic strains. The numerical predictions of the method were compared with plates molded with an experimental tool where different values for rib thickness, rib height, and base plate thickness could be set. In addition, the mold was equipped with a pressure transducer to ensure that the pressure curve was known, although simulation of the packing phase was not yet generally available. The plates were molded with a profile measuring device and good agreement with the simulation was demonstrated.

In the late 1980s and early 1990s several influential doctoral studies were carried out in the Department of Mechanical Engineering at Ohio State University. The first thesis by Ramachandra [32] involved an extensive experimental investigation of sink mark depth on a test plate with varying rib thickness and rib root radius, as well as an investigation of the effects of process parameters such as melt temperature, mold temperature, and packing pressure. The results of the study have been cited in numerous subsequent studies, and have even found their way into the calibration report of a commercial software [33]. The second thesis by Reifschneider [34] formulated a finite element approach to calculate thermal shrinkage after the pressure in the cavity had dropped to zero and the plastic was able to separate from the mold walls. This work included a comprehensive discussion of the boundary conditions of the displacement field. Finally, Beiter [35] developed a geometric index to predict the depth of sink marks for a given packing pressure. These three theses represent important contributions to the understanding of sink mark formation in injection molding.

The research was continued by Mahesh Gupta of Michigan Technological University, who used the commercial software C-Mold [36] to simulate the injection molding process for the samples produced and measured by Ramachandra [32]. The researchers transferred the temperature results calculated by C-Mold to their own two-dimensional computational domain and continued calculating until a homogeneous temperature was reached. In this way, they ensured that all details from the filling phase and flow effects were taken into account. When the pressure of the C-Mold simulation dropped to zero, a structural analysis was carried out to calculate the thermal shrinkage. The material was assumed to be elastic-perfectly plastic, with temperature-dependent elastic modulus and yield stress around the glass transition temperature. The results were compared with measurements from Ramachandra [32] and Beiter [35]. In the next article by Shi and Ramachandra [37], the study was extended to analyze the necessary size of the computational domain to obtain the same results as in the full part analysis along with the influence of the mesh size. A year later, a three-dimensional analysis was introduced in the same way [38], and the plate used in Beiter’s study [35] was analyzed with two ribs crossing in the middle.

Subsequently, several experimental studies were conducted using methods such as design-of-experiment and genetic algorithms [39,40,41,42]. New modifications of conventional injection molding techniques have been explored, including rapid heating and cooling [43], gas-assisted injection molding [44], and microcellular injection molding [45]. An interesting departure from these studies is the use of artificial intelligence trained on a large dataset of process parameters and quality assessment results [46]. Finally, in the area of numerical sink mark prediction, one study [47] considered complex elastic material properties and a cooling rate-dependent equation of state for semi-crystalline plastic, but predicted only one-third of the measured sink mark depth.

All the studies mentioned above start by calculating the deformation of the solidifying plastic when the pressure at the point of interest reaches the atmospheric value and show that such an approach can correctly predict the depth of the sink marks. However, by disregarding the exact pressure history of the filling and packing phases, the large-scale warpage is conceptually decoupled from the small-scale deformations of the parts. Moreover, all of these studies employ fairly complex elastomechanical material models such as elastic-perfectly plastic [31,36,37,38], thermo-viscoplastic [34], or complex temperature- and strain rate-dependent [47] models. In this paper, we investigate the predictive power of a method (which we refer to as the “thermoelastic method”) that uses data available in material characterization reports of a commercial injection molding simulation software, treating the solid phase as an elastic material and the molten phase as either a pure liquid or a rubbery material with an elastic modulus reduced by three orders of magnitude.

2. Thermoelastic Method

We have developed a method for predicting the local deformation of injection molded parts by dynamically calculating the displacement and stress fields of the solidifying plastic. Plastic parts are typically sheet-like, which means that some of their geometrical features can be adequately represented numerically by a two-dimensional cross-section, as already shown in Figure 1. A two-dimensional computational domain $\mathrm{\Omega }$ in the shape of the mold cavity is defined in the $xy$-plane, as shown in Figure 2. For simplicity, we assume an infinite dimension in the direction of the z-axis. We assume that plastic flow occurs only in the direction of the z-axis, so that the cross section under study is filled instantaneously at time zero. This assumption eliminates the need to model the flow. Flow-induced residual stresses are not modeled in this study. As established by Baaijens [26], they are two orders of magnitude smaller compared to thermal and pressure-induced stresses.

As soon as the cross-section under investigation is filled with hot plastic melt, heat transfer to the colder mold with homogeneous temperature begins, initiating the evolution of the temperature field of the plastic $T(x,y,t)$. It is in the nature of the injection molding process that the plastic part cools throughout the cycle, resulting in the continuous growth of a layer of solidified plastic cooled below the glass transition temperature $T_g$. Accordingly, we dynamically divide the computational domain $\mathrm{\Omega }$ into time-dependent solidified ${\mathrm{\Omega }}_s$ and molten ${\mathrm{\Omega }}_m$ subdomains with a time-dependent phase boundary $\partial {\mathrm{\Omega }}_m$ and fixed outer boundary $\partial {\mathrm{\Omega }}_s$, as shown in Figure 2. The solidified domain ${\mathrm{\Omega }}_s$ is treated as an elastic solid in all phases.

During the filling and packing phases, a hydrodynamic plastic flow in the z-direction exists in the molten ${\mathrm{\Omega }}_m$ subdomain, which is controlled externally by the volumetric filling rate or the packing pressure. The hydrodynamic pressure $p(x,y,z,t)$ corresponding to this flow is provided by the commercial software Autodesk Moldflow 2023. As verified, this pressure field can be considered homogeneous in the cross-sectional plane to a good approximation. Thus, the time-dependent pressure field $p(z,t)$ evaluated at the center of the investigated cross-section serves as a boundary condition at the phase boundary $\partial {\mathrm{\Omega }}_m$, as shown in Figure 3.

In the packing phase, the flow through a cross-section compensates for the thermal shrinkage of the downstream plastic. For a given cross-section, the packing phase ends when the hydrodynamic pressure at that cross-section drops to zero (the ambient pressure of the vented mold). This occurs when the packing flow at the set packing pressure can no longer accommodate the downstream thermal shrinkage due to the large cooling-induced viscosity increase of the upstream plastic. At this point, thermal shrinkage prevails and the pressure at this cross-section begins to drop. When it reaches the ambient value, the pressure gradient in the z-direction disappears and the flow through this cross-section is completely stopped.

In the final hydrostatic cooling phase, the pressure connection to the molding machine is interrupted and there is no longer any transport of the melt in the z-direction. The hydrodynamically isolated molten core continues to shrink thermally together with the solidified plastic, while there is no pressure that would prevent the solidified plastic from detaching from the mold walls.

In this phase, we treat the molten region ${\mathrm{\Omega }}_m$ with two approaches, either as a liquid with a shear modulus of zero as before in Figure 3, or as a rubbery material with a finite shear modulus that is three orders of magnitude lower compared to the solidified plastic, as suggested in previous studies [36,48,49] (Figure 4). In the latter case, the complete elastic problem is solved in the entire domain. In the former case, the change $\Delta p$ of the uniform pressure in the liquid core is obtained by equating the change in volume of the liquid due to thermal shrinkage and compression/dilation with the change in the volume of the cavity due to elastic displacement $\Delta \mathbf{u}$ of the solidified plastic at the $\partial {\mathrm{\Omega }}_m$ boundary and solving for $\Delta p$ by numerical iteration:

$${\int }_{{\mathrm{\Omega }}_m}dS{\alpha }_m\left(\mathbf{r}\right)\Delta T\left(\mathbf{r}\right)-\Delta p{\int }_{{\mathrm{\Omega }}_m}dS{K}_m^{-1}\left(\mathbf{r}\right)={\oint }_{\partial {\mathrm{\Omega }}_m}dl\widehat{\mathbf{n}}·\Delta \mathbf{u}(\Delta T,\Delta p).$$

(1)

In the surface integrals on the left, ${\alpha }_m$ and $K_m$, are the coefficients of thermal expansion and the bulk modulus of the liquid melt (Appendix A); $\mathbf{r}=(x,y)$. In the edge integral on the right, $\widehat{\mathbf{n}}$ is the outward normal of the ${\mathrm{\Omega }}_m$ domain and $\Delta \mathbf{u}$ corresponds to a particular field $\Delta T(x,y)$ in the solidified domain and to the sought $\Delta p$ in the liquid.

The plastic, already solidified at a given time, undergoes elastic deformation caused by the changes in the time-dependent stress field due to thermal shrinkage and/or the changes in the time-dependent pressure $p\left(t\right)$ of the molten core. As solidification progresses, the solidified and molten regions constantly change. Therefore, at each time t we are dealing with a separate elastic problem defined only during a short time interval $\Delta t$ for the region ${\mathrm{\Omega }}_s\left(t\right)$, which is solid at time t. Similarly, in the rubbery core treatment of the cooling phase we are dealing with an elastic problem defined only in the time interval $\Delta t$ for the regions ${\mathrm{\Omega }}_s\left(t\right)$ and ${\mathrm{\Omega }}_{m,r}\left(t\right)$ as they occur at time t.

The result of the elastic problem in the time interval $\Delta t$ is an increment $\Delta u_i$ of the displacement field and an increment $\Delta {\sigma }_{ij}$ of the stress field, with which the total displacement field $u_i\left(t\right)$ and the stress field ${\sigma }_{ij}\left(t\right)$ are updated:

$${\sigma }_{ij}\left(t\right)=\sum ^t\Delta {\sigma }_{ij},u_i\left(t\right)=\sum ^t\Delta u_i.$$

(2)

The stress field describes the physical state of the material, and its divergence is zero at all times. In contrast, the total displacement field relates directly to the physical strain only in the regions not traversed by the $\partial {\mathrm{\Omega }}_m$ phase boundary during evolution. Strictly speaking, in our case this is only the outer boundary. However, in the case of overmolding, for example, where cool plastic is placed in the larger mold cavity and then molded over, this “regular” region would be larger.

The crucial implication is that both the final stress state and the final shape of the cross-section generally depend on the specific time course of all relevant variables, i.e., of the fields $T(x,y,t)$ and $p\left(t\right)$, and not on the mere difference between their final and initial values, as would be the case with a constant ${\mathrm{\Omega }}_s$ and ${\mathrm{\Omega }}_{m,r}$.

The heat equation (Section 2.1) and the thermoelastic problem (Section 2.2) are solved numerically with the finite element method [50] using a suitable discretization; see Figure 5. The discretization of the elastic domains ${\mathrm{\Omega }}_s$ and ${\mathrm{\Omega }}_{m,r}$ is updated at each time step according to the change in the shape of the phase boundary $\partial {\mathrm{\Omega }}_m$. The increments of the displacement and stress fields from the individual time steps of the elastic problems are accumulated by mapping them onto a mesh of the entire domain $\mathrm{\Omega }$, which must be fine enough to capture the details of the displacement and stress field increments of all time steps. The heat equation, on the other hand, allows for a coarser mesh of $\mathrm{\Omega }$, which is used in this case.

2.1. Heat Equation

The evolution of the temperature field $T(x,y,t)$ in the entire domain $\mathrm{\Omega }$ is modeled by the heat equation

$$\rho c_p\dot{T}=\nabla ·(\lambda \nabla T),$$

(3)

where $\dot{T}\equiv \partial T/\partial t$. We assume that there is no temperature gradient along z, and consequently no heat flux in this direction. The temperature-dependent density $\rho$, specific heat $c_p$, and thermal conductivity $\lambda$ of the plastic material can be obtained from material characterization reports provided by commercial software for thousands of plastic grades. For certain materials, the thermal parameters $\lambda$ and $c_p$ are measured over a wide temperature range, as discussed in Appendix B, while for others only a single value is provided. The density is determined by the pvT relation, as described in Appendix A, and is commonly available.

A homogeneous initial condition ${T(t=0)|}_{\mathrm{\Omega }}=T_m$ is employed, where $T_m$ is the average temperature of the melt that has reached the investigated cross section. Typically, this value is a few degrees below the melt temperature, which is a process parameter. Such an average temperature result provided by commercial software takes into account the average viscous heating and conductive heat loss experienced by the plastic on its way to the investigated cross-section. In general, it would be ideal if our method had access to the entire temperature field calculated by detailed modeling of the filling phase. However, we argue that the homogeneous initial condition for the specific cross-section and process setup leads to a temperature field comparable to the Moldflow calculation (Figure 6).

Due to the imperfect contact between the polymer and the mold walls, the heat transfer through this interface must be modeled with a heat transfer coefficient h [51], which implies the boundary condition

$$\lambda {\left.\nabla T\right|}_{\partial \mathrm{\Omega }}·\widehat{\mathbf{n}}=h(T_{\mathrm{mold}}-T_{\partial \mathrm{\Omega }}),$$

(4)

where $\widehat{\mathbf{n}}$ is the outward normal of $\partial \mathrm{\Omega }$. While the appropriate value of h is widely discussed in the literature [51,52], we use the same approach as Moldflow [53] with specific values for the filling, packing, and cooling phases:

$$h_{\mathrm{filling}}=5000\frac{\mathrm{W}}{{\mathrm{m}}^2\mathrm{K}},h_{\mathrm{packing}}=2500\frac{\mathrm{W}}{{\mathrm{m}}^2\mathrm{K}},h_{\mathrm{cooling}}=1250\frac{\mathrm{W}}{{\mathrm{m}}^2\mathrm{K}}.$$

(5)

Figure 7 illustrates the evolution of the phase boundary, which initially matches the shape of the computational domain and then gradually contracts to enclose regions with larger heat capacity that take longer to cool.

2.2. Elastomechanics

The solidified/rubbery plastic is treated as an elastic material obeying Hooke’s law [54] for changes in the strain $\Delta u_{ij}(x,y,t)$ and stress $\Delta {\sigma }_{ij}(x,y,t)$ fields in individual time steps $\Delta t$ with well-defined computational domains ${\mathrm{\Omega }}_s\left(t\right)$ and ${\mathrm{\Omega }}_{m,r}\left(t\right)$:

$$\Delta {\sigma }_{ij}=K\Delta u_{kk}{\delta }_{ij}+2G\left(\Delta u_{ij}-\frac{1}{3}\Delta u_{kk}{\delta }_{ij}\right)-K\alpha \Delta T{\delta }_{ij},$$

(6)

where K and G are the bulk and shear elastic moduli, respectively, and ${\delta }_{ij}$ is the Kronecker delta. The last term represents the stress due to thermal expansion/contraction, where $\alpha$ is the coefficient of thermal expansion. The temperature-dependent $\alpha$ and K are extracted from the pvT relation (Appendix A). The shear modulus is determined from K and the Young’s modulus E, which is measured at room temperature for the solid phase and is reduced by three orders of magnitude for the rubbery material, as discussed in Appendix C.

The displacement vector $\Delta \mathbf{u}(x,y,t)$ is calculated numerically by the finite element method using a standard weak (variational) formulation. Notwithstanding the fact that the initial stress ${\sigma }_{ij}\left(t\right)$ before the onset of $\Delta \mathbf{u}$ is nonzero, the elastic problem for $\Delta \mathbf{u}$ is as usual, because the initial state at time t is an equilibrium one. This can be seen immediately from the variation of the elastic free energy density

$$\begin{array}{ll}\delta f& =\left({\sigma }_{ij}+\Delta {\sigma }_{ij}\right)\delta \Delta u_{ij}=\left({\sigma }_{ij}+\Delta {\sigma }_{ij}\right){\partial }_j\delta \Delta u_i\end{array}$$

(7)

$$\begin{array}{ll}& ={\partial }_j\left[\left({\sigma }_{ij}+\Delta {\sigma }_{ij}\right)\delta \Delta u_{ij}\right]-\left({\partial }_j{\sigma }_{ij}+{\partial }_j\Delta {\sigma }_{ij}\right)\delta \Delta u_i,\end{array}$$

(8)

where in line (7) we have used the fact that the stress tensor is symmetric. The first term of line (8) is a divergence, and is converted into a boundary contribution ${\sigma }_{ij}{\widehat{n}}_j$ (with $\widehat{\mathbf{n}}$ the outward normal) which is balanced by external forces on the boundary, as the system is in equilibrium at time t. For the same reason, in the second term of line (8), ${\partial }_j{\sigma }_{ij}=0$ everywhere. Thus, the new equilibrium condition at time $t+\Delta t$ simply requires ${\partial }_j\Delta {\sigma }_{ij}=0$, where $\Delta {\sigma }_{ij}$ is provided by Equation (6), which is the usual elastomechanical problem for $\Delta \mathbf{u}$. From the first term of line (8), we obtain the boundary condition

$$\Delta {\sigma }_{ij}{|}_{\partial }{\widehat{n}}_j=\Delta {\sigma }_{ij}^{\mathrm{ext}}{\widehat{n}}_j,$$

(9)

where $\Delta {\sigma }_{ij}^{\mathrm{ext}}$ is any change in the external stresses on the boundary within the time interval $\Delta t$.

2.2.1. Boundary Conditions

While the part is in the mold $\Delta u_z=0$ is assumed because of the shape of the part beyond the investigated cross-section that does not allow for shrinkage in the z-direction, simple ribs or a slight widening of the cross-section are sufficient to restrict the displacement along z. However, the elastic force is present in this direction, and the stress tensor has the form

$$\Delta \sigma =\left[\begin{array}{ccc}\Delta {\sigma }_{xx}& \Delta {\sigma }_{xy}& 0\\ \Delta {\sigma }_{xy}& \Delta {\sigma }_{yy}& 0\\ 0& 0& \Delta {\sigma }_{zz}\end{array}\right].$$

(10)

On the outer boundary $\partial {\mathrm{\Omega }}_s$, we use two types of boundary conditions:

• Fully constrained, meaning that while the plastic pushes against the mold, the displacement there is zero: ${\Delta \mathbf{u}|}_{\partial {\mathrm{\Omega }}_s}=0$. Consequently, the displacement is small everywhere, and most of the induced stresses cannot be relaxed. An example of the displacement and stress fields for this case is shown in Figure 8.
• Free boundary, meaning that when the contact force between the plastic and the mold reaches zero the plastic is free to detach from the mold. The boundary condition is now a zero force on the surface of the plastic, $\Delta {\sigma }_{ij}{\widehat{n}}_j{|}_{\partial {\mathrm{\Omega }}_s}=0$. The in-plane displacement allowed here relaxes most of the in-plane stress. The stress component ${\sigma }_{zz}$ cannot relax significantly because of the remaining constraint $\Delta u_z=0$. An example of the displacement and stress fields for this case is shown in Figure 9.

The actual boundary condition at the hard wall of the mold is complex, and may involve a mixture of the two cases or even a third scenario where the plastic is pressed against the mold wall but does not adhere to it, allowing tangential movement. Solidified plastic, however, usually moves away from the mold walls due to negative thermal shrinkage. Only a highly inhomogeneous thermal load could cause such global deformation that the plastic would push locally against the mold walls, as shown in Figure 9. However, limiting the local intrusion into the mold only affects the large-scale deformation of the part and not its local features such as sink marks, which are our main focus.

2.2.2. Thermal Stress Relaxability: The Role of Geometry and Melt Elasticity

In the case of the free outer boundary, much of the thermally induced stress is relaxed by the displacement. The remainder consists of residual stresses due to the inhomogeneity of the thermal load and the elastic response of the medium, which together prevent complete relaxation in the presence of general thermal inhomogeneity.

To demonstrate differences in this behavior, two cross-sections with different tendencies to displace the volume (Equation (1)) enclosed by the phase boundary during thermal shrinkage were studied:

• The circular cross-section (Figure 10a) displaces a minimal volume of solidified plastic at the phase boundary. Despite the (radial) inhomogeneities in the $\Delta T$ field, the displacement does not change the shape of the cavity, and it only shrinks according to the average effect of the temperature change. It should be noted that the coefficient of thermal expansion of a typical plastic in the molten state is about three times higher than that of the solid state; see Figure A3. Consequently, the equilibrium pressure resulting from the balance between the stiffness of the outer shell and the thermal shrinkage of the molten core (Equation (1)) is strongly negative.
• The rectangular cross-section in Figure 10b deforms in a more complex and global way. The inward bending of the thin long walls of solidified plastic displaces even more volume than the free thermal shrinkage of the molten core, resulting in a slightly positive equilibrium pressure in the core. Due to the low bending stiffness of the outer shell, a small pressure is sufficient to bring the displaced volume to an equilibrium magnitude.

The global bending deformation seen in Figure 10b is subject to global displacement of the core liquid; therefore, it depends crucially on its unobstructed flow. A finite, albeit very small, elastic modulus of the melt prevents its large-scale flow, which is required to compensate for global variations of the core cross-section. Consequently, global deformations affecting the shape of the cavity are suppressed in the rubbery case, as can be seen from the flattened bending deformation in Figure 11b. In contrast, geometries without such global deformations are not decisively affected by the weak elasticity of the melt; see Figure 11a.

Figure 12 shows the result for an even longer rectangular domain. Such a cross-section corresponds to the simplest possible injection molded part, a plate. Empirically, it is clear that the homogeneous displacement in the center, away from localized end effects, corresponds to reality.

2.3. Ejection Warpage

When a part has cooled to a uniform temperature, it no longer deforms. As long as it remains in the mold with restricted displacement in the z-direction, it has a large ${\sigma }_{zz}$ component of the stress tensor. Other stress components may be significant if the displacement of the frozen layer has been constrained by the mold walls. When the part is removed from the constraints of the mold, the boundary conditions change to those of zero forces and a new elastic equilibrium must be determined in the ejection step.

This is accomplished by solving another elastomechanical problem, analogous to the previous problems for each time step, with the boundary condition according to Equation (9) as before and by updating the displacement and stress fields according to Equation (2) as before. Thus, the boundary condition is now $\Delta {\sigma }_{ij}{{|}_{\partial }{\widehat{n}}_j=-{\sigma }_{ij}|}_{\partial }{\widehat{n}}_j$, where ${\sigma }_{ij}$ is the stress accumulated up to the event of ejection. In particular, $\Delta {\sigma }_{zz}=-{\sigma }_{zz}$, resulting in shrinkage along z. The calculated displacement field increment $\Delta \mathbf{u}$ (Figure 13) is added to the displacement field accumulated up to the ejection event to obtain the final total displacement field $\mathbf{u}$ of the plastic after ejection.

The deformation that occurs when a part is removed from the tool is commonly referred to as warpage. Commercial software usually calculates this type of deformation only once at the end of the production process. While the stress field required for this is not accumulated over the multitude of elastic problems of individual time steps, as in our thermoelastic method, it is an effective stress field constructed by various methods from simple to more sophisticated. One of the simplest is the linear shrinkage approach; this uses an isotropic stress field that, when released, causes the same local volumetric shrinkage as local cooling [18]:

$${\sigma }_{ij}(x,y)=K\frac{v\left(T(x,y)\right)-v\left(25°\mathrm{C}\right)}{v\left(T\right(x,y\left)\right)}{\delta }_{ij}.$$

(11)

The local cooling shrinkage is calculated as the relative volume change between the local temperature $T(x,y)$ at the beginning of the cooling phase and the reference room temperature $25°\mathrm{C}$. The specific volume v is provided by the equation of state (Appendix A).

The present thermoelastic approach distinguishes between deformations inside the mold, which are subject to relevant boundary conditions, and those occurring outside the mold, where free boundary conditions apply. By applying the thermal load to the developing frozen layer rather than to the fully solidified cross-section, this method takes into account local intricacies, and thereby allows the prediction of deformation patterns with small characteristic length scales.

3. Results

The final surface deformation, which is the central result of the thermoelastic method, is now demonstrated using the example of a cross-section of an experimental specimen subjected to a time-dependent pressure $p\left(t\right)$ at a packing pressure $50\mathrm{MPa}$, which is examined in detail in Section 3.4.

Figure 14 shows different deformations calculated with several variants of the method. In the region of increased thickness, a deformation in the form of a classical sink mark can be observed in all cases.

In the primary method shown in Figure 14a, the plastic is allowed to detach from the mold walls when the pressure drops.

In the simplified method shown in Figure 14b, detachment from the mold walls is restricted during all time steps. Instead of surface deformation, stress builds up, causing the surface of the part to deform only when it is removed from the mold. Although there is an equivalence between stress and displacement via Hooke’s law, this equivalence only applies to a single time step with an unchanged solidified domain. The stresses that are built up in earlier time steps when the solidified domain is thin and more deformable result in less displacement if it is not allowed at that time, but only later in a completely solidified domain. Consequently, smaller local deformation is expected when detachment is restricted, which is confirmed in Figure 14b.

In the ultimately simplified variant of the method shown in Figure 14d, the final deformation is calculated in a single step from the generally inhomogeneous isotropic stress field defined in Equation (11). This linear shrinkage approach yields a spread-out result that may predict global deformation reasonably well in the context of whole-part analysis, but fail to predict local defects. It can serve as a reference for the scope of an approach that solves only a single elastic problem.

The restricted detachment is interesting because it could be of similar accuracy to sophisticated residual stress prediction methods used by commercial software as an improvement to the isotropic shrinkage method [55]. If one wanted to limit the restricted detachment approach to a single calculation of the elastic problem, this would prevent determination of the actual instantaneous pressure in the molten core according to Equation (1). By neglecting the elastic displacement at the phase boundary, the pressure changes in the molten core are overestimated, while this is partially compensated by the underestimation of the local deformation due to the restriction of the detachment. The effect of such an approximation is demonstrated in Figure 14c. The sink mark depth is comparable to that in the primary method, but the surface opposite the sink mark deforms excessively, and the shape of the sink mark is different. Further investigation is needed to determine the validity of this approach.

3.1. Pressure of the Molten Core—Prediction of Voids

The pressure change in the molten core, $\Delta p$ in Equation (1), is integrated from the beginning of the cooling phase, and can reach considerable negative values depending on the shape of the domain (Figure 15a). A sufficiently strong negative pressure can lead to the evaporation of a component within the plastic material, creating a molding defect known as a void [56] (Figure 15b). Vaporization increases the volume of the molten core and decreases the magnitude of the pressure change, which reduces the associated sink mark. The onset of voids depends on the vapour pressure of a particular plastic grade using particular additives and the type of preprocessing, including drying, which determines the moisture content.

3.2. Showcase

To support the central claim that shrinkage, global warpage, and local deformations of a part have the same origin and can only be distinguished by characteristic length, we have calculated the influence of three different pressure traces, shown in Figure 16, on the deformation of the part cross-section in Figure 17. The original pressure trace corresponds to the $50\mathrm{MPa}$ case in Section 3.4 and serves as a reference; the other two are generated synthetically, and are generally achievable through appropriate machine settings, part design, and mold construction. The lower pressure curve has the same shape as the original, except with the values reduced by a factor of 10. The longer pressure curve is obtained by multiplying the time axis after the filling phase by $1.33$.

By separating the warpage and the sink mark in Figure 17 on the basis of their different length scales and comparing the magnitudes of all three deformation components in Figure 17a–c, it can be seen that all three components can be controlled by a single process parameter, namely, the time-dependent pressure of the molten core $p\left(t\right)$. In addition, it can be seen that the depth of the sink mark is mainly determined by the time at which the pressure drops to zero, while warpage and shrinkage are more influenced by the pressure conditions in the filling and packing phases. The empirical fact that sink marks and warpage can be controlled almost separately by the height and length of the time-dependent pressure is an additional advantage.

3.3. Difference between Liquid and Rubbery Core

Naka et al. [31] provided a profilometer measurement of the sink mark under the rib for the specific case of a $2\mathrm{m}\mathrm{m}$-thick base plate with a $2.4\mathrm{m}\mathrm{m}$-thick rib molded from acrylonitrile butadiene styrene (ABS) amorphous plastic with packing pressure of $20\mathrm{MPa}$. We replicated the part geometry in Figure 18a, estimated the feeding system from the sketches, and simulated the injection molding process with Autodesk Moldflow using the molding parameters provided by Naka to obtain the time-dependent pressure $p\left(t\right)$ of the molten core in Figure 18b, which is a necessary input for our method.

The displacement of the central cross-section was calculated by treating the molten core as a liquid and as a rubbery material. A comparison of the two approaches with Naka’s profilometer measurement is shown in Figure 19. It can be seen that the liquid approach overestimates the width of the sink mark, while the rubbery approach provides the correct depth and width of the sink mark.

If the pressure persists until the molten core is small and localized below the rib, the liquid and rubbery approaches provide similar results. However, if the pressure drops to zero early when the molten core extends over most of the part, the rubbery approach, with its finite shear modulus, prevents the thermal shrinkage load from acting on the entire molten core, instead localizing its effect and producing narrower sink marks compared to the liquid approach.

3.4. Validation

Experimental samples were molded from amorphous polycarbonate (PC) plastic Makrolon 2405 with three constant packing pressure settings of $25\mathrm{MPa}$, $50\mathrm{MPa}$, and $75\mathrm{MPa}$ maintained for $6\mathrm{s}$. The valve gates were closed before the end of packing to prevent flowback. All other parameters were set according to the material supplier’s recommendations, and the injection rate was set to achieve a $1\mathrm{s}$ filling time per $100\mathrm{m}\mathrm{m}$ flow length. The produced samples were 3D scanned, and the deformation in the direction of the surface normal was extracted on the path shown in Figure 20b.

The molding process was simulated with Autodesk Moldflow software to calculate the time-dependent pressure $p\left(t\right)$ at the investigated cross section (Figure 20a). The comparison between the surface deformation of the molded samples and its numerical prediction within the liquid core approach is shown in Figure 21. An alignment between the 3D-scanned and nominal shapes is carried out to ensure that the extracted surface deformation was symmetric and zero at the lowest point. In addition, a small parabolic displacement is added to the calculated results to match the measured global warpage, while the shape of the sink mark is left unchanged.

We can conclude that the agreement between the calculated and measured data is good. Due to the long packing phase, the molten core is localized under the rib, which makes the liquid approach suitable.

4. Discussion

In this study, we have introduced and applied a method to predict the deformation of injection molded parts. Our approach advances current research on sink mark formation by incorporating both the filling and packing phases of the injection molding process. We establish that global part warpage and localized deformations such as sink marks are due to the same cause, namely, thermal shrinkage and local pressure within the molten core.

By adopting a dual perspective of the molten core, treating it as both a liquid and a rubbery material according to previous studies [36,37,38], we demonstrate that the latter is generally essential for accurate prediction of deformations. On the other hand, the liquid approach can provide correct results if the packing phase is sufficiently long that the molten core can localize when the pressure drops to zero, i.e., at the beginning of the cooling phase.

Although the limitations of the liquid approach are evident, it remains a valuable conceptual tool. It naturally describes the buildup of negative pressure in the molten core that leads to the formation of voids. Furthermore, by studying the outcomes of this simulation model, practitioners can develop an intuitive understanding of where sink mark risks might be apparent when examining the cross-sectional shape of the molded part. This understanding is aided by the identification of stiffer and weaker regions of the solidified shell in relation to the negative pressure load of the molten core.

Looking ahead, several refinements could be incorporated into the thermoelastic method, including:

• Viscoelastic relaxation: with this modification, the value of elastic parameters could be adjusted locally based on the material history. Material parameters for such an approach are increasingly common.
• Variable heat transfer coefficient: during the cooling phase, the heat transfer coefficient could be updated at each time step based on the local distance between the plastic and the mold.
• Evaluation of local detachment: the condition for detachment from the mold walls could be evaluated locally, implementing the boundary conditions for the displacement accordingly. The numerical stability of such an approach remains an open question for further investigation.

The true utility of this method lies not in its theoretical demonstration but in its potential for practical implementation in commercial software, which would make it accessible to the plastic design community. Advances in computing power have made such an implementation feasible within existing three-dimensional simulation procedures. Warpage calculation of predefined critical areas with a limited number of discretization points could be improved with such a method. Future validation procedures would reveal the extent to which the filling and packing phases are already adequately addressed within current procedures and whether the primary benefit of this method does indeed lie in its capability to incrementally calculate the displacement of solidified plastic during the cooling phase.