**1. Introduction**

The extrudate swell flow is a well-known benchmark problem in the polymer processing field, where a free-surface and a boundary stress singularity at the die exit are present. The variety of solutions for free-surface flows is usually very limited to simple cases [1]. However, with the increasing demand for large diameters and high productivity, plastic pipes have recently been produced [2].

The polymer production process is mainly non-isothermal in nature, including heat transfer phenomena such as heating and cooling, e.g., in the injection molding [3,4] or extrusion [5–8] processes. Since the flow properties are strongly dependent on rheology and temperature, it is of great interest to understand and predict this type of heat transfer phenomenon. Although there is increasing research effort on non-isothermal free-surface flows, these studies are limited to academic purposes and problems. Extrudates in general can swell up to 15% due to thermal effects [9,10]. Phuoc and Tanner [11] developed a finite element scheme based on the Galerkin method to explore the effects of thermally induced property changes in extrusion. The swelling of self-heating extruded jets was investigated, and a new phenomenon, the so-called thermal extrudate swell, was discovered. Extruded expansions of up to 70% of the die diameter have been found in a Newtonian fluid with thermal properties similar to those of low density polyethylene. Subsequently, Karagiannis et al. [12] developed a three-dimensional (3D) non-isothermal code to study viscous

**Citation:** Fernandes, C.; Fakhari, A.; Tukovic, Ž. Non-Isothermal Free-Surface Viscous Flow of Polymer Melts in Pipe Extrusion Using an Open-Source Interface Tracking Finite Volume Method. *Polymers* **2021**, *13*, 4454. https://doi.org/10.3390/ polym13244454

Academic Editor: Alexander Malkin

Received: 17 November 2021 Accepted: 11 December 2021 Published: 19 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

free-surface flows with exponential dependence of viscosity on temperature. Again, the code was based on Galerkin's finite-element formulation, and, apart from the phenomenon of thermally induced extrudate swelling, the bending and distortion of the extrudate due to temperature differences and/or geometric asymmetries were confirmed numerically and experimentally. Finally, Spanjaards et al. [13] developed a 3D transient non-isothermal finite element code to predict the extruded shape of viscoelastic fluids emerging from an asymmetric keyhole-shaped die. A systematic study was carried out to decouple the effects of shear-thinning, elasticity, and temperature in the shape of the extrudate.

Another important effect to be considered when predicting the swelling of the extrudate is the inertia, through the dimensionless Reynolds number, *Re*, defined as the ratio of inertial forces to viscous forces within a fluid that is subjected to relative internal movement due to different fluid velocities. In Newtonian flow, the extrudate swell ratio, *χ*, defined as the ratio of the final extrudate radius to that of the die, decreases with increasing *Re* and finally achieves a contraction rather than swelling, since *χ* is less than unity [14,15]. An important class of materials commonly used in extrusion is the yield-stress fluid, which does not deform when subjected to a stress below the yield stress, τ0. Bingham [16] developed the simplest constitutive equation with a yield-stress, which was then combined with a power-law model to account for the effects of shear-thinning or shear-thickening, known as the Herschel–Bulkley model [17]. To solve fluid flows described by constitutive equations with unyielded (τ ≤ τ0) and yielded regions (τ > τ0), a special approach to overcome the numerical difficulties created by these models (the viscosity diverges to infinity as the strain rate *γ*˙ approaches zero) must be adopted. For that purpose, a common methodology is to do a regularization of the constitutive equation so that the same expression is uniformly valid for any value of *γ*˙ , i.e., to ensure that the viscosity is a continuous function of strain rate. Notice that here τ and *γ*˙ represent the magnitudes of the stress tensor τ and rate-of-strain tensor *γ*˙ , respectively. Papanastasiou [18] developed the regularization model for the Bingham constitutive equation, which was then extended to work with the Herschel–Bulkley model by Ellwood et al. [19]. A survey of regularization models used in the scientific literature can be found in Frigaard and Nouar [20]. For creeping flow conditions (*Re* 1), there are several studies where the results for the extrudate swell of materials with yield stress are presented [18,21–24]. Generally, the conclusion is that in creeping extrudate swell flow of yield-stress materials, the extrudate swell ratio decreases as yield stress increases. Additionally, it was noticed that beyond a certain value of τ0, the swell contracts slightly (*χ* < 1) and reaches a minimum, from which it starts to increase again to reach unity asymptotically. For moderate Reynolds numbers (1 < *Re* < 200), the scientific literature devoted to the study of the effects of inertia on the extrudate swell of yield-stress fluids is scarce [19,25]. In this situation, as the yield stress increases, swelling at lower Reynolds numbers and contraction at higher Reynolds numbers are reduced.

In recent decades, Computational Fluid Dynamics (CFD) simulation has expanded as a significant framework to help engineers understand and improve polymer processing techniques. Nonetheless, the current state of the art of numerical methods for complex problems, which is the case of non-isothermal free-surface extrudate swell flow of yield-stress fluids, is insufficient for practical needs. In addition, proprietary packages are expensive and cannot be modified to account for the engineer's needs. As the solver here presented is based in an open-source finite volume computational library, OpenFOAM [26], this allows more capabilities to be added to the code to take into account the physics involved in the phenomenon of interest, thereby making the simulation more realistic. Therefore, a general and freely available solver, based on the open-source framework OpenFOAM [26], for the simulation of the non-isothermal free-surface flow of yield-stress fluids would be of great importance. Hence, this work deals with the development of a solver to calculate non-isothermal free-surface flows of yield-stress fluids as an extension of the isothermal free-surface flow solver presented elsewhere [15,27]. Notice that in Fakhari et al. [15], we developed the predecessor of the current code, which was used to describe Newtonian isothermal fluid flows, and there the efficiency of this framework was demonstrated in

terms of both considering a larger computational time-step and also the smaller CPU wall time to perform each iteration of the simulation. The objective of this work is to investigate the combined effects of inertia, yield-stress, and temperature on the extrudate swell of a pipe flow. It is worth pointing out that although the melted polymer chain strongly depends on shear and stretch [28], this work focuses on the macroscale dimensional distribution of the melted polymer in the extruded swell phenomena and does not go through the microscale changes of the material such as crystallization kinetics.

The paper is organized as follows: in Section 2, we present the description of the problem, specifically, the balance equations that govern the non-isothermal free-surface flow of generalized Newtonian fluids, the boundary and initial conditions, and the Arbitrary Lagrangian–Eulerian formulation for the mesh movement. In Section 3, the numerical method used to implement the proposed model is described. In Section 4, we present the numerical results obtained for the non-isothermal extrudate swell of the Herschel–Bulkley yield-stress fluid model, taking into account the combined effects of inertia, yield-stress and temperature. Finally, in Section 5, we summarize the main conclusions of this work.

#### **2. Problem Description**

In this section, the governing equations to solve the incompressible and non-isothermal polymer-melt extrudate flow in dies of cylindrical shape for pipe production are described.

#### *2.1. Balance Equations*

The equations of motion governing the non-isothermal flow of incompressible generalized Newtonian fluids inside an arbitrary volume *V* bounded by a closed moving surface *S* include continuity, momentum, and energy equations, which can be written as:

$$\nabla \cdot \mathbf{u} = \mathbf{0},\tag{1}$$

$$\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot ((\mathbf{u} - \mathbf{u}\_{\mathcal{S}})\mathbf{u}) = -\nabla P + \nabla \cdot \mathbf{\tau},\tag{2}$$

$$
\rho \, c\_p \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) = k\_p \nabla^2 T + \mathbf{\tau} : \mathbf{D}\_\prime \tag{3}
$$

where ∇ is the Hamilton differential operator, **u** is the fluid velocity vector, **u***<sup>S</sup>* is the velocity of surface *S*, *P* is the kinematic pressure, τ is the deviatoric stress tensor, *ρ* is the fluid density, *c<sup>p</sup>* is the heat capacity, *T* is the temperature, *k<sup>p</sup>* is the isotropic thermal conductivity based on the heat flux Fourier's law, and **D** = (∇**u** + ∇**u** *T* )/2 is the strain rate tensor. Note that the last term on the right hand side of Equation (3) represents viscous heat dissipation and is written in a form that assumes all mechanical energy is dissipated as heat.

#### *2.2. Constitutive Equation*

The following constitutive equation is adopted to describe the relationship between stress and strain for the viscous fluid flow,

$$\boldsymbol{\pi} = 2\eta(\dot{\gamma})\mathbf{D}\_{\prime} \tag{4}$$

where *γ*˙ = √ 2**D** : **D** is the shear rate (invariant of the strain rate tensor **D**) and *η*(*γ*˙) is the apparent viscosity, which is described by the Herschel–Bulkley model [17] with the Papanastasiou regularization [18,19],

$$\eta(\gamma) = \min \left\{ \eta\_{0\prime} \pi\_0 \dot{\gamma}^{-1} [1 - \exp(-m\dot{\gamma})] + k \dot{\gamma}^{n-1} \right\},\tag{5}$$

where *k* is the consistency constant of proportionality, *n* is the flow index exponent, which measures the degree to which the fluid is shear-thinning or shear-thickening, *m* is a stress growth parameter, and τ<sup>0</sup> is the yield shear stress. Note that the original Papanastasiou regularization does not include the artificial upper-bounding by *η*0. However, this bounding is needed to avoid an infinite viscosity for *γ*˙ → 0 and *n* < 1. The Herschel–Bulkley model is reduced to the Newtonian model when the yield-stress is zero, i.e., τ<sup>0</sup> = 0, and the flow index exponent is one, i.e., *n* = 1.

#### *2.3. Temperature Dependency of the Apparent Viscosity*

To analyze the effect of the temperature dependency of the apparent viscosity on hydrodynamic and thermal behaviour of the flow, we keep constant the flow index *n* and the yield stress τ<sup>0</sup> and consider only a temperature-dependent consistency *k*(*T*):

$$k(T) = a\_T(T)k\_0\tag{6}$$

where *a<sup>T</sup>* denotes the shift factor and *k*<sup>0</sup> is the consistency viscosity at the reference temperature *T*0. The shift factor *a<sup>T</sup>* is defined by the Williams–Landel–Ferry (WLF) function:

$$\log(a\_T) = \frac{-c\_1(T - T\_0)}{c\_2 + T - T\_0} \,\prime \tag{7}$$

where *c*<sup>1</sup> and *c*<sup>2</sup> are material parameters.

#### *2.4. Boundary and Initial Conditions*

Considering that the fluid phases are immiscible, the fluid flow Equations (1)–(3) can be used for each phase individually, and at the interface, the proper boundary conditions must be used. First, the *kinematic condition* states that the fluid velocities on the two sides of the interface, **u**<sup>1</sup> *<sup>f</sup>* and **u**<sup>2</sup> *<sup>f</sup>* , must be continuous

$$
\mathfrak{u}\_{1f} = \mathfrak{u}\_{2f}.\tag{8}
$$

Then, from the momentum conservation law, Equation (2), follows the *dynamic condition*, which states that forces acting on the fluid at the interface are in equilibrium,

$$\left[T\_2 - T\_1\right] \cdot \mathbf{n} = \nabla\_s \sigma - \sigma \kappa \mathbf{n}\_\prime \tag{9}$$

where *T*<sup>1</sup> and *T*<sup>2</sup> are the stress tensors defined in terms of the local fluid pressure and velocity fields as *T*<sup>1</sup> = −*p*1**I** + *η*1(*γ*˙)[∇**u**<sup>1</sup> + (∇**u**1) *T* ] and *T*<sup>2</sup> = −*p*2**I** + *η*2(*γ*˙)[∇**u**<sup>2</sup> + (∇**u**2) *T* ], respectively, *σ* is the interfacial tension and ∇*<sup>s</sup>* = [**I** − **nn**] · ∇ = ∇ − **n** *∂ ∂n* is the tangential gradient operator, which appears because *σ* and **n** are defined only on the surface. From Equation (9) we derive the normal and tangential force balances [29] appropriate at a fluid–fluid interface,

$$p\_2 - p\_1 = \sigma \kappa - 2(\eta\_2(\dot{\gamma}) - \eta\_1(\dot{\gamma})) \nabla\_s \cdot \mathbf{u}\_\prime \tag{10}$$

$$\begin{aligned} \eta\_2(\dot{\gamma})[\mathbf{n} \cdot (\nabla \mathbf{u}\_l)\_2] - \eta\_1(\dot{\gamma})[\mathbf{n} \cdot (\nabla \mathbf{u}\_l)\_1] &= -\nabla\_s \sigma \\ &- \mathbf{n} (\eta\_2(\dot{\gamma}) - \eta\_1(\dot{\gamma})) \nabla\_s \cdot \mathbf{u} \\ &- (\eta\_2(\dot{\gamma}) - \eta\_1(\dot{\gamma}))(\nabla\_s \mathbf{u}) \cdot \mathbf{n}, \end{aligned} \tag{11}$$

where *κ* = −∇*<sup>s</sup>* · **n** is twice the mean curvature of the interface and **u***<sup>t</sup>* = (*I* − **nn**) · **u** is the tangential velocity component.

Dirichlet or Neumann boundary conditions for temperature *T* are specified depending on the boundary wall considered. Natural convection on boundaries are treated by setting the gradient (Neumann type boundary condition) according to

$$
\nabla T\_b = \frac{h}{k\_p} (T\_b - T\_\infty) \tag{12}
$$

where the temperature on the boundary *T<sup>b</sup>* is obtained from the previous iteration, *h* is the convection heat transfer coefficient, and *T*<sup>∞</sup> is the ambient temperature.

#### *2.5. Arbitrary Lagrange–Eulerian Formulation*

The above mathematical model, valid for arbitrary moving volume, is obtained from the corresponding material volume model using the Reynolds' transport theorem. For an arbitrary moving volume, the relationship between the rate of change of the volume *V* and the velocity **u***<sup>s</sup>* is defined by the *geometrical* (*space*) *conservation law* [30],

$$\frac{d}{dt} \int\_{V} dV - \oint\_{S} \mathbf{n} \cdot \mathbf{u}\_{s} \, dS = 0. \tag{13}$$

The problem of extrudate swelling contains moving boundaries due to the movement of the free surfaces of the extrudate. Thus, the domain is described with a mesh that is moving in time in such a way that the mesh moves with the free surfaces.

#### **3. Numerical Method**

The mathematical model for the non-isothermal free-surface flow of an incompressible generalized Newtonian fluid described in Section 2 is numerically discretized using the finite volume method [27]. First, the numerical integration in time is performed using a second-order accurate implicit method [31], referred to as the backward scheme. Next, the integral forms of the fluid flow equations are discretized in space using a second-order accurate cell-centred unstructured finite volume method. The spatial domain is discretized using a mesh, which is constituted by finite volumes (the so-called cells or elements) with an arbitrary volume *V* bounded by a closed moving surface *S*, that conserve the relevant quantities, such as mass, momentum, and energy:

$$\oint\_{S} \mathbf{n} \cdot \mathbf{u} \, dS = 0,\tag{14}$$

$$\frac{d}{dt} \int\_{V} \mathbf{u} \, dV + \oint\_{S} \mathbf{n} \cdot (\mathbf{u} - \mathbf{u}\_{S}) \mathbf{u} \, dS = - \int\_{V} \nabla P \, dV + \oint\_{S} \mathbf{n} \cdot \mathbf{\tau} \, dS,\tag{15}$$

$$\rho c\_p \left( \frac{d}{dt} \int\_V T \, dV + \oint\_S \mathbf{n} \cdot (\mathbf{u}T) \, dS \right) = \oint\_S \mathbf{n} \cdot (k\_p \nabla T) \, dS + \int\_V (\pi : \mathbf{D}) \, dV. \tag{16}$$

Detailed information of the finite volume discretization employed in the moving mesh interface tracking algorithm can be found in Tukovi´c and Jasak [27]. To summarize, the surface integrals of an integral conservation equation are transformed into sums of face integrals which together with the volume integrals are approximated to second-order accuracy by using the mid-point rule. Therefore, Equations (14)–(16) for each cell are written as,

$$\sum\_{f} \mathbf{n}\_f^n \cdot \mathbf{u}\_f^n \, \mathcal{S}\_f^n = \mathbf{0},\tag{17}$$

$$\frac{3\mathbf{u}\_P^n V\_P^n - 4\mathbf{u}\_P^o V\_P^o + \mathbf{u}\_P^{oo} V\_P^{oo}}{2\Delta t} + \sum\_f (\dot{m}\_f^n - \dot{\mathcal{U}}\_f^n) \mathbf{u}\_f^n = -(\nabla P)\_P^n V\_P^n + \sum\_f \mathbf{n}\_f^n \cdot \mathbf{\tau}\_f^n S\_{f'}^n \tag{18}$$

$$\begin{aligned} \rho\_P(c\_p)\_P \frac{3T\_P^\eta V\_P^\eta - 4T\_P^\eta V\_P^\rho + T\_P^{\eta \nu} V\_P^{\rho \alpha}}{2\Delta t} + \sum\_f \rho\_f(c\_p)\_f \mathbf{n}\_f^\eta \cdot \mathbf{u}\_f^\eta T\_f^\eta S\_f^\eta &= \sum\_f (k\_p)\_f \mathbf{n}\_f^\eta \cdot (\nabla T)\_f^\eta S\_f^\eta \\ &+ (\mathbf{\tau}\_P^\eta : \mathbf{D}\_P^\mu) V\_P^\eta \end{aligned} \tag{19}$$

where ∆*t* is the time-step, the subscripts *P* and *f* represent the cell-center and face-center values at cell with volume *VP*, the superscripts *n*, *o* and *oo* represent values evaluated at the new time instance *t <sup>n</sup>* and two previous time instances *t <sup>o</sup>* and *t oo* = *t <sup>o</sup>* <sup>−</sup> <sup>∆</sup>*t*. Finally, the cell-face mass flux *m*˙ *n <sup>f</sup>* = **n** *n f* · **u** *n f S n <sup>f</sup>* must satisfy the discretized mass conservation law (Equation (17)), and the face volume flux *U*˙ *<sup>n</sup> <sup>f</sup>* must satisfy the discretized geometrical conservation law (Equation (13)). Equation (17) is linearized for a pressure correction

according to the consistent non-iterative PISO algorithm [32–34]. The algorithm can be summarized as follows [35]:


The Poisson-type equation for pressure is solved with a conjugate gradient method with Cholesky preconditioner and the velocity and temperature linear systems are solved using BiCGstab with an Incomplete Lower-Upper (ILU) preconditioning [38–40].

#### **4. Results and Discussion**

In this section, we present the validation and assessment of the newly developed moving mesh finite-volume interface tracking solver that is able to efficiently handle inelastic non-Newtonian matrix-based free-surface flows with non-isothermal effects. In this work, the fluid rheology is described by the Herschel–Bulkley constitutive equation with the Papanastasiou regularization. The newly developed solver is tested against fluid

flow simulations in an axisymmetric domain geometry typical of pipe extrusion. The first case study is devoted to the extrudate swell of a Bingham fluid, and subsequently, the extrudate swell of Herschel–Bulkley shear-thinning and shear-thickening flows is presented. These studies aim to verify the solver's capabilities to accurately predict the extrudate swell of inelastic non-Newtonian matrix-based fluids. Additionally, the effects of inertia and yield stress on the extrudate swell are investigated. Finally, the non-isothermal extrudate swell effects are studied for a Bingham fluid at moderate Reynolds number (*Re* = 10), in order to test the robustness of the newly developed numerical algorithm, specifically for non-isothermal calculations.

#### *4.1. Problem Domain and Meshes*

The benchmark case study that will be discussed is the axisymmetric extrudate swell of non-Newtonian inelastic fluids. A schematic representation of the computational flow domain, the boundary faces, and the discretization mesh for the initial time step (*t* = 0) and at a steady state is shown in Figure 1. The level of mesh refinement used in the numerical studies carried out in this study corresponds to the same level employed in mesh M5 of our previous work [15], which resulted from a mesh convergence analysis. Polar coordinates are employed for the description of the axisymmetric flow domain, thus **x** = (*r*, *z*). The half width of the axisymmetric channel is denoted as *R*, which is considered to be the scaling length. The inlet plane is taken sufficiently far upstream from the exit so that the flow is fully developed with a mean velocity *U*. In the axisymmetric domain, the two lateral boundary sides are considered to be wedge patches (i.e., the cylinder is specified as a wedge of small angle, e.g., 5◦ , and a thickness cell running along the plane of symmetry, encompassing one of the coordinate planes). At the bottom, the axis of symmetry is considered as empty patch. At the solid die wall, the no-slip (tangential velocity is zero) and no-penetration (normal velocity is zero) conditions are imposed for velocity, zero-gradient pressure, and a fixed value temperature. At the free-surface, the kinematic condition, Equation (8), and the dynamic condition, Equations (10) and (11), are imposed for pressure and velocity, along with natural convection, Equation (12), for temperature. Finally, the outflow plane is taken sufficiently far downstream of the die exit so that the flow is uniform. The die exit of the axisymmetric domain is located at *x* = 5*R* from the inlet, and the outflow is located at *x* = 25*R* from the die exit.

The dimensionless numbers governing the flow are the Reynolds number,

$$\text{Re} = \frac{\rho l I^{2-n} R^n}{k} \,\text{,}\tag{20}$$

which is the ratio of inertial forces to viscous forces within a fluid that is subject to relative internal movement due to different fluid velocities; the Bingham number,

$$Bn = \frac{\pi\_0 R^n}{k! I^n},\tag{21}$$

which is the ratio of yield stress to viscous stress and describes the extent to which the controllable yield stress can exceed the viscous stress; the dimensionless growth exponent,

$$M = \frac{m\mathcal{U}}{\mathcal{R}},$$

that, by using a regularized constitutive equation such as the Papanastasiou model, it will determine up to which convergent results can be obtained (in our simulations *M* ≥ 500 for 0 < *Bn* ≤ 10); and the Prandtl number,

$$Pr = \frac{kc\_p}{k\_p} \left(\frac{\mathcal{U}}{\mathcal{R}}\right)^{n-1} \text{ \textit{\textsuperscript{n-1}}}\tag{23}$$

which is the ratio of moment diffusivity (kinematic viscosity) and thermal diffusivity of a fluid, expressing the relationship between the movement quantity diffusion and the heat quantity diffusion within the fluid itself.

**Figure 1.** Schematic representation of (**a**) the axisymmetric extrudate swell domain geometry and boundary faces (**b**) of an indicative discretization mesh at the initial time step *t* = 0, and (**c**) at steady state.

#### *4.2. The Effects of Inertia and Yield-Stress in the Extrudate Swell of Bingham Fluids*

Initially, the combination of inertia and yield-stress effects on the extrudate swell ratio obtained for the isothermal flow of Bingham fluids (i.e., *n* = 1) is studied. The swell ratio is defined as the height of the free-surface away from the die exit, where the plug flow has been established, divided by the die radius, i.e., *χ* = *h*0/*R*. As shown in Figure 2, the extrudate swell ratio obtained by the newly developed solver is very close to the results presented by Kountouriotis et al. [41] at *Re* = 1, 5, and 10. For the lowest Reynolds number (*Re* = 1), as the yield stress effect is enhanced (i.e., increase of the *Bn* number), the extrudate swell shrinks steeply for *Bn* > 0.1, and above a critical value of *Bn* > 5, the swell contracts and becomes smaller than unity. At *Re* = 5 the extrudate swell ratio increases with increasing *Bn*, reaching a maximum at *Bn* ≈ 1. Subsequently, the extrudate swell ratio decreases and becomes lower than unity for *Bn* > 5. At *Re* = 10, the extrudate swell ratio is below one for the lowest *Bn* numbers (*Bn* < 0.5), and then it enlarges, being higher than one as *Bn* enhances, reaching a maximum at *Bn* ≈ 2. After that the swell ratio *χ* diminishes and again becomes lower than unity for *Bn* > 5. Regardless of the *Re* number used in the simulations, for *Bn* ≥ 5, the extrudate swell ratio is very similar on each curve, which means that the yield stress effects are predominant compared to the inertial effects.

**Figure 2.** Steady-state extrudate swell ratio *χ* for the isothermal extrudate swell of Bingham fluids (*<sup>n</sup>* <sup>=</sup> 1) at *Re* <sup>=</sup> {1, 5, 10} and *Bn* <sup>=</sup> {10−<sup>3</sup> , 10−<sup>2</sup> , 10−<sup>1</sup> , 0.5, 1, 5, 10}. Solid lines represent the results obtained by Kountouriotis et al. [41], and the symbols represent the results obtained by the newly developed interface tracking code.

In Figure 3, the magnitude of the polymer velocity vector is shown for *Re* = 1, 5, and 10 and *Bn* = 0.001 and 10. At *Bn* = 0.001 (negligible yield-stress effects), the maximum value of the magnitude of the polymer velocity vector is twice the inlet velocity and the enhancement of inertia effects through the increase in the *Re* number leads to the reduction of the die swell ratio, mimicking the behavior of the extrudate swell for a Newtonian fluid, as presented by Fakhari et al. [15]. When the highest *Bn* number is considered (*Bn* = 10), the maximum value of the magnitude of the polymer velocity vector is 1.3 times the inlet velocity. This means that the increase in the yield stress promotes the retardation of the flow, changing from parabolic shape profile to plug-flow. The consequence is that the extrudate swell ratio decreases and, in fact, contracts. Note that the contour plots of the magnitude of the polymer velocity vector are very similar at this higher Bingham number, where increasing inertia effects do not change the shape of the flow field.

Figure 4 shows the contours of the magnitude of the polymer stress tensor at steadystate for *Re* = 1, 5, and 10 at *Bn* = 0.001 and 10. When the yield stress effects are negligible (*Bn* = 0.001), the maximum value of the stress tensor magnitude occurs for the lowest Reynolds number (*Re* = 1), corresponding to the highest extrudate swell ratio. Increasing the inertia effects decreases the localized maximum value of the magnitude of the polymer stress tensor in the top right corner of the die wall, leading to a decrease in the extrudate swell ratio. On the other hand, when the yield stress effects dominate the flow (i.e., *Bn* = 10), the maximum value of the magnitude of the polymer stress tensor is constant for all *Re* numbers from 1 to 10, being three orders of magnitude lower than the case with *Bn* = 0.001.

**Figure 3.** Contours of the magnitude of the polymer velocity vector at steady-state for the isothermal extrudate swell of Bingham fluids (*n* = 1) when *Re* = 1 (**top**), *Re* = 5 (**middle**) and *Re* = 10 (**bottom**), and *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**).

**Figure 4.** Contours of the magnitude of the polymer stress tensor at steady-state for the isothermal extrudate swell of Bingham fluids (*n* = 1) when *Re* = 1 (**top**), *Re* = 5 (**middle**) and *Re* = 10 (**bottom**), and *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**).

#### *4.3. The Effects of Inertia and Yield-Stress in the Extrudate Swell of Herschel–Bulkley Fluids*

In this section, the effects of inertia and yield-stress are investigated on the isothermal extrudate swell of Herschel–Bulkley flows (*n* 6= 1). For this purpose, the simulations were carried out at the same *Re* and *Bn* numbers as those presented in Section 4.2, but with a flow index exponent representative of shear-thinning, *n* = 0.5, and shear-thickening, *n* = 1.5, behaviors. As shown in Figure 5, the extrudate swell ratio obtained by the newly developed solver is again very close to the results presented by Kountouriotis et al. [41] for both *n* = 0.5 and *n* = 1.5 at *Re* = 1, 5, and 10. The general trend of the shear-thinning extrudate swell ratio (see Figure 5a) at the different *Re* numbers is similar to the behavior shown for the Bingham flow (*n* = 1) in Figure 2. In Figure 5a, we can see that, in general, for *n* = 0.5, the extrudate swell ratio is lower than the ones obtained for *n* = 1 and *n* = 1.5. This behavior allows us to conclude that the shear-thinning nature of the Herschel–Bulkley fluid helps to reduce the swelling of the polymer. On the other hand, for the shear-thickening behavior (*n* = 1.5), the extrudate swell ratio at different *Re* and *Bn* numbers is the highest, as shown in Figure 5b. Additionally, the variation in the extrudate swell ratio with the increase in the yield stress effects is more abrupt for this shear-thickening fluid, reaching a maximum of 20% at *Re* = 1.

**Figure 5.** Steady-state extrudate swell ratio *χ* for the isothermal extrudate swell of Herschel–Bulkley fluids at *Re* <sup>=</sup> {1, 5, 10} and *Bn* <sup>=</sup> {10−<sup>3</sup> , 10−<sup>2</sup> , 10−<sup>1</sup> , 0.5, 1, 5, 10} with (**a**) *n* = 0.5 and (**b**) *n* = 1.5. Solid lines represent the results obtained by Kountouriotis et al. [41], and the symbols represent the results obtained by the newly developed interface tracking code.

Figure 6 shows the contours of the magnitude of the polymer velocity vector for the isothermal extrudate swell of shear-thinning Herschel-Bulkley fluids (*n* = 0.5), when *Re* = 1, 5 and 10 and *Bn* = 0.001*s* and 10. At *Bn* = 0.001 (negligible yield stress effects), the maximum value of the magnitude of the polymer velocity vector is approximately 1.7 times the inlet velocity for 1 ≤ *Re* ≤ 10, which is 15% smaller than the value obtained for the Bingham flow (*n* = 1), leading to variations of 7% in the extrudate swell ratio when increasing *Re* from 1 to 10. On the other hand, when the highest *Bn* number is considered (*Bn* = 10), the maximum value for the magnitude of the polymer velocity vector is 1.1 times the inlet velocity, which is again 15% smaller than the value obtained for the Bingham flow (*n* = 1). However, for this strongly dominated yield stress flow, the extrudate swell ratio is kept constant for 1 ≤ *Re* ≤ 10 and again with a plug-flow velocity profile leading to swell contraction (*χ* < 1).

**Figure 6.** Contours of the magnitude of the polymer velocity vector at steady-state for the isothermal extrudate swell of shear-thinning Herschel–Bulkley fluids (*n* = 0.5) when *Re* = 1 (**top**), *Re* = 5 (**middle**) and *Re* = 10 (**bottom**), and *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**).

In Figure 7, we show the contour plots of the magnitude of the polymer stress tensor at steady-state for the isothermal extrudate swell of shear-thinning Herschel–Bulkley fluids (*n* = 0.5) when *Re* = 1, 5, and 10 and *Bn* = 0.001 and 10. When the yield stress effects are negligible (*Bn* = 0.001), the maximum value of the magnitude of the polymer stress tensor decreases with the increase of inertia from *Re* = 1 to 10, while at the highest *Bn* number, the inertia do not have influence on the maximum value of the magnitude of the polymer stress tensor. In general, for both of the *Bn* numbers, the magnitudes of the polymer stress tensor are approximately one third of the ones obtained for the Bingham fluid (*n* = 1).

Figure 8 displays the contours of the magnitude of the polymer velocity vector for the isothermal extrudate swell of shear-thickening fluids (*n* = 1.5), when *Re* = 1, 5, and 10 and *Bn* = 0.001 and 10. As can be seen by the velocity contours at *Bn* = 0.001, the maximum value for the magnitude of the polymer velocity vector is 2.2 times the inlet velocity for 1 ≤ *Re* ≤ 10, which is 10% larger than the value obtained for the Bingham flow (*n* = 1), leading to variations of 17% in the extrudate swell ratio when *Re* is increased from 1 to 10. At the highest Bingham number (*Bn* = 10), increasing the inertia effects does not change the extrudate swell ratio and the maximum value of the polymer velocity vector magnitude. The later remains 1.5 times the inlet velocity, which is approximately 13% larger than the one obtained for the Bingham fluid (*n* = 1).

**Figure 7.** Contours of the magnitude of the polymer stress tensor at steady state for the isothermal extrudate swell of shear-thinning Herschel–Bulkley fluids (*n* = 0.5) when *Re* = 1 (**top**), *Re* = 5 (**middle**) and *Re* = 10 (**bottom**), and *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**).

**Figure 8.** Contours of the magnitude of the polymer velocity vector at steady state for the isothermal extrudate swell of shear-thickening Herschel–Bulkley fluids (*n* = 1.5) when *Re* = 1 (**top**), *Re* = 5 (**middle**) and *Re* = 10 (**bottom**), and *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**).

Finally, in Figure 9, we show the steady-state contour plots of the magnitude of the polymer stress tensor for the isothermal extrudate swell of shear-thickening Herschel– Bulkley fluids (*n* = 1.5) when *Re* = 1, 5, and 10 and *Bn* = 0.001 and 10. As before, the stress contours follow the same trend as in the other fluids, i.e., the maximum magnitude of the polymer stress tensor decreases with an increase in inertia at *Bn* = 0.001, while it is constant for the yield-stress-dominated flow (*Bn* = 10), regardless of the enhancement of inertia from *Re* = 1 to 10. Note also that the maximum magnitude of the polymer stress tensor for the shear-thickening Herschel–Bulkley fluid is around twice larger than the value obtained for the Bingham fluid at both *Bn* numbers (*Bn* = 0.001 and 10).

**Figure 9.** Contours of the magnitude of the polymer stress tensor at steady-state for the isothermal extrudate swell of shear-thickening Herschel–Bulkley fluids (*n* = 1.5) when *Re* = 1 (**top**), *Re* = 5 (**middle**) and *Re* = 10 (**bottom**), and *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**).

#### *4.4. Non-Isothermal and Yield Stress Effects in the Extrudate Swell of Bingham Fluids*

In this section, we study the effects of temperature and yield stress in the extrudate swell ratio of Bingham fluids (*n* = 1) at Reynolds number *Re* = 10. We consider two different scenarios for the die wall temperature *Tw*, one where *T<sup>w</sup>* < *Tinlet* (cold wall) and another where *T<sup>w</sup>* > *Tinlet* (hot wall). We will then examine the thermally induced swelling behavior as the Bingham number increases for these two configurations.

The thermal and physical properties of a typical polystyrene [7] used in the simulations are listed in Table 1. An important issue for modeling the cooling of the polymer when it is extruded is the definition of the boundary condition at the polymer and air interface. In this work, we employed a defined convective heat flux as given in Equation (12). Additionally, the influence of the temperature on the rheological behavior of the material is controlled by the WLF equation for the shift factor employed in the temperature-dependent consistency

parameter *k* (see Equations (6) and (7)). Typical extreme sets of WLF parameters (*c*1, *c*2) are (4.54, 150.36), leading to thermorheological coupling [42].

**Table 1.** General conditions used in the non-isothermal simulations of the extrudate swell of Bingham fluids (*n* = 1).


In this case study, the ratio of moment diffusivity and thermal diffusivity of the fluid, defined by the dimensionless Prandtl number, is given by *Pr* = 0.7 (see Equation (23)). A small *Pr* number means that heat diffuses very easily compared to moment.

In Figure 10, we show the extrudate swell ratio obtained by the newly developed solver for the non-isothermal flow of Bingham fluids at *Re* = 10 and 0.001 ≤ *Bn* ≤ 10, with a cold (blue symbols) and hot (red symbols) die wall. For the cold wall, the extrudate swell ratio has a decreasing monotonic behavior, but being always larger than unity, meaning that the polymer expands after the die wall. Note that the decreasing monotonic behavior is related to the enhancement of the yield stress effects. Specifically, the extrudate swell ratio *χ* slightly decreases from *Bn* = 0.001 to *Bn* = 0.1, and then it falls with a sharper slope, having its smallest value around 1.06 at *Bn* = 10. In contrast with the cold die wall, the extrudate swell ratio for the heating case study follows the same trend as the isothermal studies, where we first see an increase in the extrudate swell ratio followed up by a decrease due to the larger yield stress effects. Only at *Bn* = 2 and 3, the polymer expands after the die wall, but for all the other *Bn* numbers, the polymer melt contracts, which is a similar behavior to the shear-thinning isothermal extrudate swell at *Re* = 10 presented in Section 4.3. Note that in Labsi et al. [43], it is shown that neglecting the viscosity's temperature dependency leads to undervaluing hydrodynamic properties, especially in the cooling case. If we compare the results shown in Figures 2 and 10 at *n* = 1 and *Re* = 10, we see that the extrudate swell ratio for the case neglecting viscosity's temperature dependency (Figure 2) is smaller than the one obtained when viscosity's temperature dependence is taken into account (Figure 10), and this behavior is more noticeable for the cooling case, corroborating the conclusion presented in Labsi et al. [43].

Figure 11 shows the magnitude of the polymer velocity vector at *Re* = 10 for *Bn* = 0.001 and 10, and for the cold and hot die wall case studies. For *Bn* = 0.001 (negligible yield-stress effects), the maximum value of the magnitude of the polymer velocity vector is 30% higher in the cold die wall case study and 5% lower in the hot die wall case study, than the maximum values obtained for isothermal conditions. When the higher *Bn* number is considered (*Bn* = 10), the maximum value of the magnitude of the polymer velocity vector is 46% higher in the cold die wall case study and 8% lower in the hot die wall case study, than the maximum values obtained for isothermal conditions. Notice also that for the cold die wall case study at *Bn* = 10 the velocity profile has changed from a plug-flow shape in isothermal conditions to a parabolic profile shape now. This means that the non-isothermal effects are stronger than the yield stress effects for the cold die walls case study at *Re* = 10. For the hot die wall case study at *Bn* = 10, the velocity profile maintains the plug-flow shape profile, leading to a reduction in the extrudate die-swell ratio.

**Figure 10.** Steady-state extrudate swell ratio *χ* for the non-isothermal axisymmetric extrudate swell of Bingham fluids (*<sup>n</sup>* <sup>=</sup> 1) at *Re* <sup>=</sup> 10 and *Bn* <sup>=</sup> {10−<sup>3</sup> , 10−<sup>2</sup> , 10−<sup>1</sup> , 0.5, 1, 2, 3, 4, 5, 10}.

**Figure 11.** Contours of the magnitude of the polymer velocity vector at steady-state for the nonisothermal asymmetric extrudate swell flow of Bingham fluids (*n* = 1) at *Re* = 10 for *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**), when *T<sup>w</sup>* < *Tinlet* (**top**) and *T<sup>w</sup>* > *Tinlet* (**bottom**).

In Figure 12, we show contours of the steady-state dimensionless temperature field at *Re* = 10 for *Bn* = 0.001 and 10, and for the cold and hot die wall case studies. As can be seen in Figure 12, the temperature variation ∆*θ* for the case where cold die walls are

used in the simulations is approximately twice the one obtained when hot die walls are employed, which will not induce a high swelling ratio for the latter case (hot die walls). This is due to the fact that the material flows toward the cooler (more viscous) region as the flow rearranges in the extrudate, which will give higher swelling ratio of the cooler side (as shown in the cold die walls).

**Figure 12.** Contours of dimensionless temperature at steady-state for the non-isothermal asymmetric extrudate swell flow of Bingham fluids (*n* = 1) at *Re* = 10 for *Bn* = 0.001 (**left**) and *Bn* = 10 (**right**), when *T<sup>w</sup>* < *Tinlet* (**top**) and *T<sup>w</sup>* > *Tinlet* (**bottom**).

#### **5. Conclusions**

A numerical algorithm able to solve non-isothermal and inelastic non-Newtonian free-surface flows based on the Arbitrary Lagrangian–Eulerian (ALE) formulation was presented and implemented using the finite-volume method. The implementation was performed in the open-source OpenFOAM framework [26], where the interface is tracked in a semi-implicit manner, which allowed robust and stable deformations of the interface.

The newly developed algorithm was assessed in terms of accuracy for isothermal flow simulations of the axisymmetric extrudate swell of Bingham and Herschel–Bulkley fluids. The effects of inertia and yield stress on the extrudate swell ratio computed by the newly developed algorithm were studied and compared with results found in the scientific literature for the range of the dimensionless Reynolds (*Re*) and Bingham (*Bn*) numbers as follows: 1 ≤ *Re* ≤ 10 and 0.001 ≤ *Bn* ≤ 10. Additionally, the contours for the magnitude of the polymer velocity vector and stress tensor fields are also shown and discussed. For the isothermal Bingham flows, the extrudate swell ratio was found to vary by approximately 13% from the lowest to the highest *Re* number when the yield stress effects are negligible (*Bn* = 0.001). Additionally, at higher *Bn* numbers, the yield stress effects are dominant and the extrudate swell ratio is equal for all the *Re* numbers tested, being lesser than unity, which means that the swell contracts. For the isothermal Herschel– Bulkley flows of the shear-thinning fluid, the extrudate swell ratio was found to vary by approximately 8% from the lowest to the highest *Re* number when the yield stress effects are negligible (*Bn* = 0.001). On the other hand, at higher *Bn* numbers, the yield stress effects are dominant and the extrudate swell ratio is equal for all the *Re* numbers tested, being also lesser than unit as in the Bingham flow case, meaning that the swell contracts. In the case of using a shear-thickening fluid, the extrudate swell ratio was found to vary by approximately 17% from the lowest to the highest *Re* number when the yield stress effects are negligible (*Bn* = 0.001). Additionally, at higher *Bn* numbers, the extrudate swell

ratio approaches the unitary value, which means that the polymer melt neither contracts nor expands.

Finally, the newly developed algorithm was applied to study the non-isothermal flow of axisymmetric extrudate swell using Bingham fluids at *Re* = 1. The effects of temperature and yield stress on the extrudate swell ratio computed by the newly developed algorithm were analyzed for the range of the dimensionless Bingham (*Bn*) numbers, 0.001 ≤ *Bn* ≤ 10, and two different configurations, one representing a cold die wall and the other an hot die wall. Additionally, the contours for the magnitude of the polymer velocity vector and temperature fields are shown and discussed. For the cold die wall it was found that the extrudate swell ratio has a monotonic decreasing behavior but always with a value greater than unity for 1 ≤ *Bn* ≤ 10. However, for the hot die wall, the extrudate swell ratio first increases for 0.001 ≤ *Bn* ≤ 2 and then decreases for 2 ≤ *Bn* ≤ 10, being greater than unity only at *Bn* = 2 and 3.

In summary, the results presented here show that the newly developed interface tracking code can accurately predict the non-isothermal extrudate swell of inelastic non-Newtonian matrix-based fluids. The code that was implemented here is being currently extended to handle viscoelastic non-Newtonian fluid flow calculations.

**Author Contributions:** Conceptualization, C.F.; Formal analysis, C.F.; Investigation, C.F. and A.F.; Methodology, C.F. and A.F.; Resources, C.F.; Software, C.F. and Ž.T.; Supervision, C.F. and Ž.T.; Validation, C.F. and A.F.; Writing—original draft, C.F. and A.F.; Writing—review & editing, C.F., A.F., and Ž.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to acknowledge the funding by FEDER funds through the COM-PETE 2020 Programme and National Funds through FCT—Portuguese Foundation for Science and Technology under the projects UIDB/05256/2020 and UIDP/05256/2020.

**Acknowledgments:** The authors also acknowledge the support of the computational clusters Search-ON2 (NORTE-07-0162-FEDER-000086) the HPC infrastructure of UMinho under NSRF through ERDF and FCT I.P. through the Advanced Computing Project CPCA/A00/6057/2020 using the Minho Advanced Computing Center (MACC).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Chao-Ming Lin \* and Yun-Ju Chen**

Department of Mechanical and Energy Engineering, National Chiayi University, Chiayi 600355, Taiwan; yunju.chen.zoe@gmail.com

**\*** Correspondence: cmlin@mail.ncyu.edu.tw

**Abstract:** Plastic is an attractive material for the fabrication of tubular optical instruments due to its light weight, high strength, and ease of processing. However, for plastic components fabricated using the injection molding technique, roundness and concentricity remain an important concern. For example, in the case of a telecentric lens, concentricity errors of the lens barrel result in optical aberrations due to the deviation of the light path, while roundness errors cause radial stress due to the mismatch of the lens geometry during assembly. Accordingly, the present study applies the Taguchi design methodology to determine the optimal injection molding parameters which simultaneously minimize both the overall roundness and the overall concentricity of the optical barrel. The results show that the geometrical errors of the optical barrel are determined mainly by the melt temperature, the packing pressure, and the cooling time. The results also show that the optimal processing parameters reduce the average volume shrinkage rate (from 4.409% to 3.465%) and the average deformations from (0.592 mm to 0.469 mm) of the optical barrel, and the corresponding standard deviation values are reduced from 1.528% to 1.297% and from 0.263 mm to 0.211 mm, respectively. In addition, the overall roundness and overall concentricity of the barrel in the four planes are positively correlated.

**Keywords:** plastic optical barrel; injection molding; roundness; concentricity; Taguchi method

#### **1. Introduction**

Plastic injection molding is a fast and economical process for fabricating optical components with high precision, excellent performance, and good strength-to-weight properties. However, the quality of injection molded parts is highly dependent on the choice of processing parameters, including the melt temperature, mold temperature, filling time, packing time, packing pressure, injection pressure, and so on. As a result, a proper control of the processing conditions is essential [1–7].

Compared to conventional lenses, in which the magnification varies with the distance between the lens and the object, telecentric lenses have a constant field of view at all distances from the lens. As a result, they eliminate the parallax error inherent in traditional fixed-focal length lenses and, therefore, find widespread use in machine vision-based systems where precise and repeatable measurements are required, such as metrology, microlithography, semiconductor manufacturing, and so on [8–12]. Figure 1 presents a simple schematic illustration of a typical coaxial bilateral telecentric optical system consisting of a telecentric barrel, a coaxial light source, a holder, and an optical imaging system.

Among these components, the telecentric barrel plays a critical role in ensuring that the light emitted by the coaxial light source follows the preset path. In particular, the roundness and concentricity of the barrel must be strictly controlled in order to ensure a proper placement and alignment of the internal lens [13–18]. The lens barrel is generally fabricated from plastic material due to the latter's light weight and ease of manufacturing. However, as shown in Figure 1, the barrel has a varying tube diameter and an asymmetric geometry (due to the presence of the holder). Thus, controlling the injection molding

**Citation:** Lin, C.-M.; Chen, Y.-J. Taguchi Optimization of Roundness and Concentricity of a Plastic Injection Molded Barrel of a Telecentric Lens. *Polymers* **2021**, *13*, 3419. https://doi.org/10.3390/ polym13193419

Academic Editors: Célio Bruno Pinto Fernandes, Alexandre M. Afonso, Luís L. Ferrás and Salah Aldin Faroughi

Received: 21 September 2021 Accepted: 2 October 2021 Published: 5 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

parameters in such a way as to minimize the roundness and concentricity errors caused by the shrinkage and deformation of the barrel during the molding process represents a significant challenge [19,20].

**Figure 1.** Coaxial bilateral telecentric optical system.

In general, the Taguchi method and mold flow analysis provide a convenient and cost-effective approach for optimizing the processing conditions employed in the injection molding process. Accordingly, the present study employs a hybrid approach consisting of the Taguchi design method and mold flow simulations to determine the optimal settings of the main injection molding parameters (i.e., the injection pressure, the packing pressure, the melt temperature, the mold temperature, and the cooling time) for simultaneously minimizing both the roundness and the concentricity errors of the telecentric barrel. Having determined the optimal processing conditions, a further investigation is performed to examine the correlation between the overall roundness and the overall concentricity of the barrel and the effects of the optimal processing conditions on the average volume shrinkage rate of the barrel following removal from the mold [21–25].

#### **2. Theoretical Analysis**

In the present study, the mold flow analysis is performed using Moldex3D computer aided engineering (CAE) simulations. The simulations assume a contact interface between the part's surface and the mold wall and separate the warpage analysis process into two parts, namely, in-mold deformation during the packing and cooling stages and free deformation following ejection from the mold. For the geometric accuracy requirements of the coaxial bilateral telecentric barrel, the final deformation analysis of the cured part from the temperature after demolding to the room temperature will be calculated based on the analysis method of the roundness and concentricity. The related theoretical calculations are explained as follows.

#### *2.1. Flow Analysis during Filling Stage*

Using the Moldex3D solid simulation, the polymer melt flow develops during the filling stage of the injection molding process where the melt flow is assumed to be incompressible. The polymer melt is assumed to be Generalized Newtonian Fluid (GNF). Therefore, the non-isothermal 3D flow motion can be mathematically described by the mass, momentum conservation, and energy conservation equations, which can be written as follows [26]:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \rho u = 0 \tag{1}$$

$$\frac{\partial(\rho u)}{\partial t} + \nabla \cdot (\rho uu - \sigma) = \rho \text{g} \tag{2}$$

$$
\rho c\_p \left( \frac{\partial T}{\partial t} + \mu \cdot \nabla T \right) = \nabla \cdot (k \nabla T) + \eta \dot{\gamma}^2 \tag{3}
$$

where *u* is velocity, *p* is the pressure, *ρ* is the density, *c<sup>p</sup>* is the heat capacity, *η* is the viscosity, . *γ* is the shear rate, *k* is the heat conductivity, and *σ* is the stress tensor; this can be expressed as follows:

$$\boldsymbol{\sigma} = -p\mathbf{I} + \eta \left(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^{\mathrm{T}}\right) \tag{4}$$

Considering the constitutive equation for the general polymer materials, the Modified-Cross viscosity [27] model with Arrhenius temperature is used to describe the rheological property of the polymer melt.

$$\eta\left(T,\dot{\gamma}\right) = \frac{\eta\_0(T)}{1 + \left(\eta\_0(\dot{T})\gamma/\tau^\*\right)^{1-n}}\tag{5}$$

$$
\eta\_0(T) = B e^{\left(\frac{T\_k}{T}\right)}.\tag{6}
$$

where *η* is the viscosity, *η*<sup>0</sup> is the melt viscosity under zero-shear-rate conditions, *τ* ∗ is the parameter that describes the transition region between zero shear rate and the power law region of the viscosity curve, *n* is the Power Law index, and B is the consistency index. The Modified-Cross viscosity model includes the Newton's fluid interval and the power-law shear thinning interval. When the shear rate approaches zero, this model predicts the zero-shear rate viscosity *η*0; when the shear rate is large, it predicts the power-law behavior. The *τ* ∗ in this model is a constant, which physically represents the critical shear stress value for the transition from Newtonian fluid to power law fluid. Compared with the other model, the Modified-Cross model requires fewer parameters and can capture the dependence of viscosity on the shear rate. Therefore, the Modified-Cross model is often used in commercial simulation software, just like the Moldex3D software in this study.

A volume fraction function, *f*, is introduced to track the advance of the melt front. Here, *f* = 0 is defined as the air phase and *f* = 1 as the polymer melt phase. Hence, the melt front is located within cells with 0 < *f* < 1. The advancement of f over time can be expressed as the following transport equation:

$$\frac{\partial f}{\partial t} + \nabla \cdot (u f) = 0 \tag{7}$$

After the part is ejected from the mold, a free thermal shrinkage happens due to the temperature and pressure difference. The warpage analysis assumes the mechanical properties are elastic. The stress–strain equilibrium equations enable us to solve the problems.

#### *2.2. Shrinkage and Warpage*

The temperature and pressure changes which occur during the injection molding process result in corresponding changes in the specific volume and density of the polymer. These changes lead in turn to a warpage of the molded component as it cools from the melt condition to the solid condition. The part additionally undergoes volume shrinkage during the molding process and following its removal from the mold. During the packing stage, the shrinkage reduces as the packing pressure and packing time increase.

For semi-crystalline polymers, the shrinkage behavior mainly depends on the degree of crystallization. If the mold temperature is low and the cooling rate is high, it is not easy to crystallize, but there is a small shrinkage; on the other hand, if the mold temperature is high and the cooling rate is low, the macromolecular chain has enough relaxation time and easily to form crystals. The amount of shrinkage will naturally increase.

For isotropic materials, the linear shrinkage is one-third the volumetric shrinkage (see Equation (8) below). However, in the injection molding process, the orientation effect of the polymer forming and the shrinkage behavior are both constrained by the mold wall. As a result, the shrinkage deformation exhibits a non-isotropic behavior and the linear shrinkage in the part thickness is thus governed by Equation (9) [28–33].

$$\mathbf{S}\_L = \mathbf{1} - \left(\mathbf{1} - \mathbf{S}\_V\right)^{1/3} \approx \left(\frac{1}{3}\right) \mathbf{S}\_V \tag{8}$$

$$S\_L \approx 0.9 - 0.95S\_V \tag{9}$$

where *S<sup>L</sup>* is the linear shrinkage rate, and *S<sup>V</sup>* is the volume shrinkage rate.

#### *2.3. Stress Analysis after the Demolding Stage*

In this stage, the plastic forming part is no longer restricted to the mold after demolding and is in the free shrinking stage. The free volume shrinkage of the molded component following its removal from the mold depends mainly on the thermal stress induced by the difference between the temperature of the part and that of the environment. If the shrinkage stress exceeds the mechanical strength of the part, the part undergoes warpage. Conversely, if the plastic part is sufficiently strong to resist the thermal stress, the part retains its original geometry and dimensions. However, shrinkage voids may still be formed within the plastic component, which degrade the mechanical properties of the part and may lead to cracks and breakage under the effects of an external force. In the warpage analysis of the Moldex3D solid model, the assumptions are as follows: the material property is linear and elastic; there is a small amount of strain; the behavior is approximately steady; and the plastic part is elastically deformed. The governing equations for the material behavior in the warpage analysis can thus be expressed as follows [34]:

$$
\sigma\_{\dot{l}\dot{l},\dot{l}} + f\_{\dot{l}} = 0 \tag{10}
$$

$$
\sigma\_{\rm ij} = \mathbb{C}\_{\rm ijkl} \left( \varepsilon\_{kl} - \varepsilon\_{kl}^{0} - \mathfrak{a}\_{kl} \cdot \Delta T \right) + \sigma\_{\rm ij}^{F} \tag{11}
$$

$$
\varepsilon\_{i\bar{j}} = \left( u\_{i,\bar{j}} + u\_{j,\bar{i}} \right) / 2 \tag{12}
$$

where *σij* is the stress tensor, *σ F ij* is the initial stress induced by the flow, *εij* is the infinitesimal elastic strain, *ε* 0 *ij* is the initial strain from the P-*v*-T relationship, *Cijkl* is the elastic material stiffness, *αkl* is the coefficient of linear thermal expansion, and ∆*T* is the temperature difference.

#### *2.4. Roundness Evaluation*

The most common methods for determining roundness errors include the Least Squares Circle (LSC) method, the Minimum Zone Tolerance Circle (MZC) method, the Maximum Inscribed Circle (MIC) method, and the Minimum Circumscribed Circle (MCC) method [35]. Figure 2 illustrates the LSC method, in which the center of the circle is first determined by identifying the circular contour which minimizes the sum of squared error (SSE) between the inner and outer radii of the interior surface (shown in red in Figure 2). This center point is then used to draw the circumscribed and inscribed circles of the barrel interior surface, respectively (see two black lines in Figure 2). Finally, the roundness of the circle (∆*Zq*) is quantified as the radial distance (*Rmax*–*Rmin*) between them [36]:

$$\text{Roundness} = \Delta Zq = R\_{\text{max}} - R\_{\text{min}} \tag{13}$$

**Figure 2.** Least squares circle method for determination of roundness (∆*Zq*) at specified Z-plane.

In practice, the center point of the least square circle is unique, and its accuracy depends on the number of measurement points [37–40]. The overall roundness including n planes can be defined as

$$\text{Roundness}\_{\text{overall}} = \sum\_{i=1}^{n} \left[ (\Delta Zq)\_i \right]^2 / n \tag{14}$$

After the above calculation, one can get the center of each contour (*Xc, Yc*), the roundness of each contour, and the overall roundness of all contours after the injection molding processing.

#### *2.5. Concentricity Evaluation*

Concentricity refers to the deviation of the center of a circle or center of a cylinder from the center of the reference form. It is generally evaluated as either the axis concentricity tolerance (with the tolerance zone centered on the axis of the reference form) or the point concentricity tolerance (with the reference point taken as the center of the circle). For either method, the concentricity measurement process involves establishing the coordinates of the required checking plane, measuring the position of the center of the contour circle that needs to be compared after setting the datum, and calculating the distance between the original center of the circle and the center of the actual contour. As with the roundness evaluation described above, the accuracy of the concentricity tolerance process also increases with an increasing number of measurement points. Figure 3 illustrates the concentricity evaluation process for the case where the reference coordinates (*X,Y*) are set as (0,0) and the coordinates of the fitted circular contour are denoted as (*Xc, Yc*). The concentricity *d* at specified Z-plane is then evaluated simply as [41].

$$\text{Concentricity} = d = \sqrt{(X\_{\text{c}} - X)^2 + (Y\_{\text{C}} - Y)^2} \tag{15}$$

**Figure 3.** Calculation of concentricity (d) at specified Z-plane.

The overall concentricity including n contours can be defined as

$$\text{Concentricity}\_{\text{Overall}} = \sum\_{i=1}^{n} \left[ d\_i \right]^2 / n \tag{16}$$

#### **3. Methods and Procedures**

The geometry model of the telecentric lens barrel and positions of the tracked nodes during the simulations were defined in **Rhinoceros**. The model was then imported into **Moldex3D** to design the mold cavity and the gate, runner and cooling system, as well as to perform the molding flow simulations. The simulations considered the use of **PA66** polymer material with the properties shown in Table 1 as the feedstock material. The total warpage was obtained directly from the output results of the mold flow analysis, while the roundness was calculated based on the distance between each tracked node and the offset center, and the concentricity was computed as the shortest distance between the offset circle center and the original axis.


**Table 1.** Material properties of PA66 (TECHNYL A 216, Solvay Engineering Plastics; Source: Moldex3D material library).

**Figure 4.** (**a**) Viscosity and (**b**) P-*v*-T properties of PA66 material (Source: Moldex3D material library).

#### *3.1. Material Characteristics and P-v-T curves of PA66*

Figure 4a shows the relationship between the viscosity of PA66 and the shear rate at different temperatures as the basis for subsequent mold flow analysis, and the Modified-Cross viscosity [27] model (See Equations (5) and (6)) with Arrhenius temperature is used to describe the rheological property of the polymer melt. Figure 4b shows the P-*v*-T relationship diagram of PA66. In the process of plastic processing, the plastic undergoes a very rapid cooling process under the temperature and pressure controlled by the molding process and changes from a molten state to a solid state. Usually, the volume changes greatly, and a simple comparison is no longer possible. To describe the capacity constant, the relationship between specific volume/pressure/temperature characteristics (P-*v*-T) is

determined to calculate the degree of compression of the material in the packing stage, as well as the shrinkage rate and shrinkage warpage of the final plastic part after ejection.

The Modified Tait Model II [26] is used to describe the P-*v*-T relationship of semicrystalline materials (PA66) and is also the recommended P-*v*-T model in Moldex3D.

$$v(T, P) = v\_0(T) \left[ 1 - \text{Cln} \left( 1 + \frac{P}{B(T)} \right) \right] + v\_l(T, P) \tag{17}$$

where *v*(*T*, *P*) is the specific volume; *v*<sup>0</sup> is the specific volume at zero gauge pressure; *T* is the temperature; *P* is the pressure; and *C* is the constant 0.0894.

$$\vartheta\_0(T) = \begin{cases} \begin{array}{c} b\_{1S} + b\_{2S}\overline{T} \\ b\_{1L} + b\_{2L}\overline{T} \end{array} , \text{if } T \le T\_{\text{fl}} \\ \text{if } T > T\_{\text{fl}} \end{cases} \tag{18}$$

$$B = \begin{cases} \begin{array}{l} b\_{3S} \exp\left(-b\_{4S}\overline{T}\right) \\ b\_{3L} \exp\left(-b\_{4L}\overline{T}\right) \end{array} , if \ T \le T\_{\text{f}}\\ \text{ (19)}\\ B = \begin{array}{l} \text{ (10)}\\ \text{ (10)} \end{array} , if \ T > T\_{\text{f}} \end{cases} \tag{19}$$

$$v\_l(T, P) = \begin{cases} \ \ \ b\_l \exp\left(b\_8 \overline{T} - b\_9 P\right) & \text{, if } T \le T\_l \\\ \ \ 0 & \text{, if } T > T\_l \end{cases} \tag{20}$$

$$
\overline{T} = T - b\_5 \tag{21}
$$

$$T\_{\!\!\!\!=} = b\_{\!\!\!=} + b\_{\!\!\!=}P \tag{22}$$

where *v<sup>t</sup>* is the value for semi-crystalline resins only applies to temperatures below the transition temperature; T<sup>t</sup> is used to characterize the abrupt viscosity change of the material around its transition temperature; 13 parameters (*b*1*S*, *b*2*S*, *b*3*S*, *b*4*S*, *b*1*L*, *b*2*L*, *b*3*L*, *b*4*L*, *b*5, *b*6, *b*7, *b*8, *b*9) are data-fitted coefficients. With only linear P-*v*-T transitions, b7, b<sup>8</sup> and b<sup>9</sup> are for amorphous materials.

#### *3.2. Modeling of Analyzed Product*

The simulations considered a coaxial telecentric lens barrel with the dimensions and geometry shown in Figure 5.

**Figure 5.** Designed geometry of telecentric lens (unit: mm).

Figure 6a,b show the gating and cooling system models used in the simulations. As shown in Figure 6a, the gating system was designed with four runners to accommodate the large component size and thin wall thickness (3 mm). Moreover, the cooling water runner system was numerically designed to fit snugly around the outside surface of the plastic barrel, and a baffle-type water runner was used to prevent internal heat accumulation (see Figure 6b).

**Figure 6.** (**a**) Gating and (**b**) cooling systems for injection molding process.

#### *3.3. Taguchi Design Method*

Figure 7 presents a flowchart of the hybrid Taguchi/CAE optimization process performed in the present study to identify the plastic injection molding processing parameters which minimize the overall roundness and overall concentricity of the optical barrel. As shown, the process commenced by constructing the numerical model described in the previous section and establishing the build surface and solid mesh. Having chosen suitable signal-to-noise (S/N) ratios for evaluating the quality of each simulation outcome, the Taguchi design processes were defined by establishing the control factors and level settings of interest. For some specific quality requirements such as deformation, warpage, shrinkage, weld line, air trap, roundness, concentricity, etc., each of these quality characteristics will have different influential processing factors. Therefore, in order to find the best combination of parameters, the Taguchi method is usually used to screen the most influential factors. This method is to utilize the statistical operation of the orthogonal array (OA) to find the optimal parameter combination. In Taguchi method, OA is a general partial factorial design. It is based on an orthogonal design matrix, allowing users to consider selected subsets of multi-factor combinations at multiple levels. Orthogonal arrays are balanced to ensure that all levels of all factors are considered equally in statistics. Other less

influential parameters adopt the recommended values of polymer materials or injection molding machines. Generally speaking, the processing factors that have an influence on deformation are related to temperature and packing. After preliminary evaluation and calculation, five control factors were chosen, namely, (A) the injection pressure, (B) the packing pressure, (C) the melt temperature, (D) the mold temperature, and (E) the cooling time. As shown in Table 2, each of the five control factors was assigned four different level settings.

**Figure 7.** Flowchart showing main steps in Taguchi/CAE optimization process.


**Table 2.** Control factors and level settings used in Taguchi simulations.

Thus, the Taguchi simulations were configured in an L16(4<sup>5</sup> ) Orthogonal Array (OA), as shown in Table 3.

**Table 3.** Taguchi analysis results for overall roundness and overall concentricity.


In the present study, the aim of the optimization process was to minimize the overall roundness and overall concentricity of the selected planes in the plastic barrel. Hence, in evaluating the quality of the solutions obtained from each simulation run in the OA, the smaller-the-better S/N ratio was adopted for both quality measures, i.e.,

$$S/N = -10\log\left(\frac{1}{n}\sum\_{i=1}^{n} y\_i^2\right) \tag{23}$$

where *y<sup>i</sup>* is the roundness or concentricity and *n* is the number of measured points in the simulation trial.

#### *3.4. Least Squares Circle Method for Evaluation of Roundness and Concentricity*

The roundness and concentricity computations were performed at four planes distributed along the barrel length, namely, Z = 0, Z = 57.75, Z = 82.3, and Z = 117.7 (mm), respectively, as shown in Figure 8. In determining the roundness using the LSC method (see Section 2.4), the center of the least squared error circle was determined using the function [42]:

$$f(\mathbf{x}, \mathbf{y}) = \min \sum\_{i=1}^{n} \left[ r(\mathbf{x}, \mathbf{y}) - \mathbb{R} \right]^2 \tag{24}$$

where *r*(*x*, *y*) is the distance between the measured point and the known center of the circle, (*x*, *y*) are the coordinates of the measured point, *n* is the number of measured points, and *R* is the radius of the least square circle [43]. For each run in the OA array, the displacements of the measurement nodes (see Figure 8) were obtained and used to obtain the center point (Xc, Yc, Zc) and radius Rc (See Table 4) of the corresponding least square circle. As described in Section 2.4, the roundness is denoted by ∆*Zq*.

**Figure 8.** Measurement nodes used for roundness and concentricity evaluation at different Z-planes.

#### **4. Results and Discussion**

After Taguchi's optimization calculation, including roundness, concentricity and correlation analysis, the relevant analysis results will be confirmed and discussed. Table 3 shows the average S/N values of the overall roundness and overall concentricity obtained in each of the 16 runs in the OA over the four measurement planes (Z1~Z4). The table also shows the S/N values obtained under the standard injection molding conditions for the injection machine and molding material (as prescribed by the manufacturer). Finally, the table shows the S/N values for the overall roundness and overall concentricity obtained under the optimal settings of the five control factors (see Section 4.1 below).

### *4.1. Factor Rank Analysis and Optimal Process Parameters*

Figure 9a,b show the Taguchi response graphs for the overall roundness and overall concentricity, respectively. Referring to Figure 9a, it is seen that the optimal overall roundness is obtained using factor level settings of A3, B4, C1, D3, and E4, i.e., an injection pressure of 220 MPa, a packing pressure of 240 MPa, a melt temperature of 275 ◦C, a mold temperature of 90 ◦C, and a cooling time of 17 s. Furthermore, the simulation results show that the overall roundness is dominated by the melt temperature (Rank 1), packing pressure (Rank 2), and cooling time (Rank 3) in sequence. By contrast, the injection pressure and mold temperature, with smaller S/N ranges of 0.020909 dB and 0.010688 dB, respectively, have only a relatively minor effect on the overall roundness. As shown in Table 3, the S/N value of the overall roundness under the optimal processing conditions (26.755 dB) is 0.283 dB higher than that of the barrel produced under the standard processing conditions (26.472 dB). Moreover, the S/N value is also higher than that produced in any of the simulation runs in the OA. In other words, the effectiveness of the optimized parameter design in minimizing the overall roundness of the molded plastic barrel is confirmed.

Figure 9b shows the Taguchi response graph for the overall concentricity of the molded barrel. It is seen that the optimal overall concentricity is again obtained using control factor level settings of A3, B4, C1, D3, and E4. The overall concentricity is determined mainly by the packing pressure (Rank 1), melt temperature (Rank 2), and cooling time (Rank 3). The injection pressure and mold temperature once again have only a minor effect on the overall concentricity. Referring to Table 3, it can be seen that the optimal processing conditions increase the S/N ratio (19.331 dB) by 1.053 dB compared with that obtained under the standard processing conditions (18.278 dB). In addition, the S/N ratio of the optimized design is higher than that obtained in any of the 16 runs of the OA. Thus, the effectiveness of the optimal processing conditions in improving the overall concentricity of the barrel is confirmed. Notably, the results presented in Table 3 show that the optimal values of the overall roundness and overall concentricity, respectively, are obtained using the same control factor level settings. In other words, the optimal design enables the simultaneous optimization of both the overall roundness and the overall concentricity.

The results presented in Figure 9a show that the packing pressure, melt temperature, and cooling time have similar S/N values, i.e., 0.146, 0.177, and 0.129 dB, respectively. In other words, all three factors exert a similar effect on the overall roundness of the molded barrel. However, the overall concentricity is dominated by a major factor, namely, the packing pressure (S/N = 0.997 dB) (see Figure 9b) and two moderate influence factors, namely, the melt temperature (S/N = 0.412 dB) and cooling time (S/N = 0.222 dB).

For the influence of packing pressure, when the mold cavity is completely filled with plastic melt, and the plastic melt will change from high temperature and high pressure to low temperature and low pressure. Due to changes in the temperature and pressure of the plastic melt, the final filled part may have obvious shrinkage in the mold cavity. Therefore, in order to overcome the shrinkage problem, the plastic melt in the runner will be continuously filled into the mold cavity when the filling stage is completed. This is called the packing stage of injection molding. In the packing stage, the inside of the mold cavity will reach the highest pressure, and the plastic melt will continue to solidify where it contacts the lower temperature mold wall. The packing process should continue

until the injection gate is solidified. Generally speaking, increasing the packing pressure or extending the packing time will delay the curing time of the plastic melt, which will promote the dispersion of the pressure in the plastic part and reduce the volume shrinkage. Excessive packing pressure is likely to cause factors such as difficulty in demolding, high residual stress, burrs and flash. On the contrary, insufficient packing pressure will lead to larger volume shrinkage and voids and other defects.

**Figure 9.** S/N response for (**a**) Overall roundness and (**b**) Overall concentricity.

The effects of the melt temperature and cooling time mainly affect the geometric deformation of the injected part. Deformation is the most important factor that simultaneously affects the optimization of both roundness and concentricity. Generally, the control factors affecting plastic deformation are temperature and cooling. The melt temperature can determine the difference between the surface compressive stress and the internal tensile stress of the injected part during the curing process. In addition, the cooling time is mainly

because the semi-crystalline polymer needs sufficient time to crystallize during the cooling process to reduce the residual stress and the shrinkage.

#### *4.2. Correlation between Overall Roundness and Overall Concentricity*

Comparing the results presented in Figure 9a,b, it can be seen that the overall roundness and overall concentricity have identical trends in terms of their dependency on the level settings of each control factor. Furthermore, for both properties, the packing pressure, melt temperature and cooling time exert the greatest effect on the simulation outcome, while the injection pressure and mold temperature have only a minor effect. Figure 10 shows the results obtained when plotting the overall roundness values in Table 3 against the corresponding overall concentricity values. Applying a regression analysis technique to the simulation data, the correlation coefficient is determined to be R<sup>2</sup> = 0.7159 (Or R = 0.846). Considering the general correlation evaluation, the correlation coefficient of the two variables is greater than 0.7, which can be regarded as highly correlated. In other words, the overall roundness and overall concentricity are positively related to one another, which explains why they respond in a similar manner to changes in the injection molding conditions and can be simultaneously optimized using the same control factor level settings.

**Figure 10.** Analysis diagram of the correlation between overall roundness and overall concentricity.

#### *4.3. Deformation and Shrinkage of Plastic Barrel*

As shown in Figure 1, a holder structure is attached to the side of the lens barrel in order to support the barrel during use and maintain the coaxial condition of the light as it passes through the barrel. However, the addition of the holder induces a deformation of the molded barrel since the greater thickness of the holder structure relative to that of the barrel results in a corresponding reduction in the local cooling rate. Observing the left-hand schematics in Figures 11 and 12, which show the barrel produced under

the standard processing conditions, the local reduction in the cooling rate results in the formation of two regions of high-volume shrinkage due to the difference in cooling rates of the outer and inner regions of the holder structure, respectively. The greater volume shrinkage rate (Figure 11) then causes the narrower portion of the barrel to deform in the direction of the holder structure (Figure 12). However, as shown in the right-hand schematics in the two figures, the optimal processing conditions suppress the local volume shrinkage effect and reduce the barrel deformation accordingly. Figure 13 compares the overall volume shrinkage rates of the barrels produced using the standard and optimized processing conditions, respectively. (Note that in the ideal case, the volume shrinkage is equal to zero.) A detailed inspection shows that the optimal processing conditions reduce the overall average shrinkage rate and standard deviation from 4.409% to 3.465% and 1.528% to 1.297%, respectively. Similarly, Figure 14 compares the overall deformations of the barrels produced using the standard and optimized processing conditions and shows that the optimal processing conditions reduce the average deformations and standard deviation from 0.592 mm to 0.469 mm and from 0.263 mm to 0.211 mm, respectively.

**Figure 11.** Cross-sectional shrinkage of final part processed using standard conditions (**left**) and optimal conditions (**right**). (Volume shrinkage: %).

#### *4.4. Manufacturing and Processing Implications of Present Results*

Table 4 shows the eccentric coordinates, least square circle radius, roundness, and concentricity values of the ideal telecentric barrel (original design) and barrels produced under the standard and optimal processing conditions, respectively.

**Figure 12.** Side-view deformation of final part processed using standard conditions (**left**) and optimal conditions (**right**). (Displacement enlarged by factor of 10 for visualization purposes).

**Figure 13.** Shrinkage improvement of optimal process compared to standard process.

**Figure 14.** Deformation improvement of the optimal process compared to the standard process.



Interestingly, the concentricity of the barrel produced under standard processing conditions at planes Z<sup>2</sup> and Z<sup>3</sup> (see Figure 8) is better than that of the barrel produced under the optimal process conditions at the same planes. It is speculated that this may be due to the holder. However, the overall roundness of the barrel produced using the optimal processing conditions is better than that of the barrel produced using the standard processing conditions at all of the measurement planes.

In general, the present results show that to improve the overall roundness and overall concentricity of the optical telecentric barrel, it is necessary to reduce the thermally induced residual stress to the greatest extent possible. The Taguchi optimization results suggest that this can best be achieved using an appropriate injection speed, increasing the packing pressure, extending the cooling time, reducing the product thickness difference, using an appropriate material temperature, reducing the mold temperature to improve the difference with the room temperature, and appropriately selecting the gate design.

#### **5. Conclusions**

In this study, the coaxial telecentric lens barrel was analyzed by considering the tolerances of the roundness and concentricity. Mold flow technology combined with the Taguchi design method were introduced to explore the roundness and concentricity arising from the shrinkage and deformation in the injection molding process. Five control factors and four levels of orthogonal array tables were selected for the Taguchi analysis. The results show that an appropriate selection of processing factors and levels can effectively optimize the roundness and concentricity of the lens barrel injection molding process. The simulation results support the following main conclusions.


**Author Contributions:** Conceptualization, C.-M.L., Y.-J.C.; formal analysis, C.-M.L., Y.-J.C.; writing original draft preparation, C.-M.L., Y.-J.C.; data curation, Y.-J.C.; writing—review and editing, C.-M.L.; supervision, C.-M.L.; funding acquisition, C.-M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science and Technology of Taiwan, ROC, under Grant Numbers MOST 109-2221-E-415-001-MY3 and MOST 109-2221-E-415-002-MY3.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors gratefully acknowledge the financial support provided to this study by the Ministry of Science and Technology of Taiwan, ROC.

**Conflicts of Interest:** The authors declare no conflict of interest.
