A CFD-driven surrogate-based machine learning framework for multi-criteria quality optimization of family-mold processes

Abstract

This study presents a machine learning-driven framework for the constrained multi-criteria optimization of family-mold processes, where the simultaneous production of geometrically dissimilar parts creates complex flow imbalances. Dimensional deformation (warpage) and shear stress are identified as critical quality characteristics affecting the functional performance and structural integrity of molded products, while product weight is treated as a practical manufacturing constraint to ensure material efficiency. Preliminary experimental trials are conducted to verify process feasibility and to determine realistic operating ranges for key injection molding parameters. Based on these ranges, a face-centered central composite design (FCCCD) is constructed, and high-fidelity computational fluid dynamics (CFD) data are generated using Moldex3D analysis. A Kriging-based surrogate model is then developed to represent the complex and nonlinear relationships between process parameters and quality responses. The surrogate model is integrated with the Non-dominated Sorting Genetic Algorithm II (NSGA-II) to perform constrained multi-criteria optimization, where warpage and shear stress are simultaneously minimized under a prescribed weight constraint. The optimization process yields Pareto-optimal process parameter combinations that effectively balance dimensional accuracy and stress reduction while satisfying material usage requirements. Selected Pareto-optimal solutions are validated through detailed CFD simulations, yielding average relative errors of $1.56\mathrm{\%}$ and $3.58\mathrm{\%}$ for warpage and shear stress, respectively, demonstrating the predictive accuracy and robustness of the proposed framework. Overall, this study establishes an effective data-driven optimization strategy for quality-oriented design and decision support in intelligent family-mold processes.

Introduction

Injection molding has established itself as a cornerstone manufacturing process in modern industry owing to its ability to produce large volumes of plastic components with high dimensional accuracy, excellent surface finish, and consistent quality, while accommodating complex geometries and tight tolerances.^1,2 The injection molding process typically consists of four sequential stages, namely plasticizing, filling, packing/holding, and cooling, each of which plays a critical role in determining the final part quality. ³ Moreover, the high level of automation and repeatability inherent to injection molding enables stable mass production with reduced labor cost and material waste, which has led to its widespread adoption across automotive, electronics, medical devices, and consumer product manufacturing sectors.

Leveraging the inherent advantages of injection molding, mold configurations can generally be categorized into single-cavity and multi-cavity systems, where the latter is widely adopted to improve productivity in mass production. Conventional multi-cavity molds are designed to manufacture identical parts simultaneously, whereas family molds represent a special class of multi-cavity molds that enable the concurrent production of multiple parts with different geometries or functional requirements within a single molding cycle.^4–6 By consolidating multiple products into one tooling system, family-mold technology can significantly reduce tooling investment, machine setup time, and operational cost, while enhancing equipment utilization and shortening production lead time, making it particularly attractive for small- to medium-batch production and multi-product manufacturing environments.

The quality of injection-molded products, including mechanical strength, dimensional accuracy, and residual stress, is influenced not only by material properties but also by product geometry and processing conditions. However, since the geometry is typically constrained by design requirements and difficult to modify, adjusting process parameters offers a more flexible and practical approach to controlling part quality and manufacturing performance. ⁷ Critical molding parameters such as melt temperature, mold temperature, filling time, and packing pressure play a decisive role in determining product quality while simultaneously affecting productivity, cycle time, and energy efficiency. In industrial practice, these parameters are often determined based on engineering experience and subsequently tuned through trial-and-error, an approach that is both time-consuming and unreliable due to the complex thermo-viscoelastic behavior of polymer melts.^8,9

To mitigate the reliance on such costly trial-and-error cycles, high-fidelity Computational Fluid Dynamics (CFD) simulations have emerged as an indispensable tool. In the broader landscape of modern engineering research, CFD simulations serve as a fundamental cornerstone for characterizing complex fluid-thermal interactions across diverse manufacturing domains.^10–12 Specifically in the realm of polymer processing, the implementation of CFD analysis facilitates a detailed characterization of highly non-linear physical interactions and thermo-mechanical behaviors. However, the high computational cost associated with these simulations presents a new challenge; a trial-and-error approach within the simulation environment itself remains highly inefficient and time-consuming when navigating a high-dimensional design space. Without a structured search strategy, the extensive resources required for each simulation run can become a significant bottleneck, particularly in complex family-mold systems where inter-cavity interactions and inherent flow imbalances must be precisely mapped. Consequently, the limitations of both experience-driven and exhaustive simulation methods have motivated the development of systematic optimization strategies to efficiently identify optimal processing conditions and enhance both product quality and manufacturing efficiency.

Traditionally, the optimization of injection molding process parameters has been predominantly conducted using classical statistical and experimental approaches. For example, Banh et al. ¹³ applied a Grey-based Taguchi method to optimize process parameters in family-mold processes, focusing on reducing sink marks and shrinkage. Their results showed that classical statistical optimization techniques can effectively improve part quality and consistency in family mold applications. Hentati et al. ¹⁴ combined the Taguchi approach with computer-aided engineering simulation to investigate the optimization of injection molding process parameters for polycarbonate/acrylonitrile-butadiene-styrene parts. Similarly, Berihun and Bogale ¹⁵ employed a Taguchi method integrated with grey relational analysis to analyze the influence of injection molding parameters on shrinkage characteristics of PET plastic preforms, reporting a noticeable reduction in shrinkage under optimized conditions. In addition, Kumar and Bairwa ¹⁶ used the Taguchi experimental design to support the optimization of injection molding parameters related to defect reduction in polypropylene containers. The analysis suggested the dominant role of melt temperature compared with other processing variables. Furthermore, Phan et al. ¹⁷ utilized a Taguchi L16 orthogonal array integrated with multi-objective grey relational analysis to minimize weld-line defects and enhance dimensional accuracy in polycarbonate electronic housings. Overall, the reviewed studies suggest that traditional statistical optimization techniques remain applicable for systematically analyzing injection molding processes and identifying influential parameters.

More recently, data-driven optimization strategies based on artificial intelligence and machine learning have been increasingly adopted to address the growing complexity of injection molding processes. By learning nonlinear relationships from limited data and coupling predictive models with advanced optimization algorithms, these approaches offer improved modeling flexibility and computational efficiency compared to traditional techniques. For instance, El Ghadoui et al. ¹⁸ combined a backpropagation neural network and a genetic algorithm for global optimization of process parameters in plastic injection molding. Their study addressed multiple quality- and efficiency-related objectives, including cycle time, material usage, and energy consumption. Similarly, Xu et al. ⁷ proposed a hybrid artificial neural network and particle swarm optimization framework to optimize injection molding process parameters for improved impact mechanical performance. Their approach models the nonlinear relationship between process parameters and mechanical response to reduce maximum von Mises stress. In another study, Liu et al. ¹⁹ presented an integrated optimization approach for stress in injection-molded parts by considering molding, assembly, and service stages. A Kriging surrogate model combined with the expected improvement criterion was employed to optimize the maximum service stress. Along a different line of research, Gim et al. ²⁰ introduced an in-mold condition-centered optimization framework for injection molding based on explainable artificial intelligence. Their method employs machine learning to model the relationships between process parameters, in-mold conditions, and part quality, enabling interpretable and robust process optimization. In addition, Zhao et al. ⁹ investigated an optimization strategy for injection molding weight control by reformulating the process parameter optimization problem as a classification task. A support vector classifier was used to identify the decision boundary associated with the standard product weight, and particle swarm optimization was applied to tune the classifier hyperparameters. The results indicated that the method can achieve satisfactory accuracy with a relatively small number of experimental samples.

Although extensive studies have investigated process parameter optimization in injection molding using surrogate models and evolutionary algorithms, a significant portion of the literature has primarily focused on single-cavity or multi-cavity molds with geometrically identical cavities. While specialized research has begun to address family mold systems, such as the investigation by Huang et al. ⁶ which provides valuable insights into assembly behavior through a single-objective lens, such approaches may not fully reflect the competitive trade-offs inherent in multi-product layouts. Although multi-criteria optimization is widely applied in injection molding, its application frequently involves simplified configurations that may offer limited consideration for the intricate inter-cavity interactions found in family molding. In these systems, dissimilar geometries result in highly coupled flow, packing, and shrinkage behaviors, which add substantial complexity to quality control. Consequently, developing a multi-criteria optimization framework specifically tailored to the unique coupling of family-mold processes provides a valuable advancement for the field.

To address this research gap, the present study develops a machine learning-driven digital twin framework that integrates high-fidelity CFD simulations, surrogate modeling, and evolutionary optimization to explicitly capture complex cavity interactions. The proposed approach enables Pareto-based optimization of warpage and shear stress under a material-efficiency-driven weight constraint, thereby providing an intelligent decision-support tool for quality-oriented family mold design in smart manufacturing environments.

Methodology

Problem statement

This study addresses the multi-criteria optimization problem associated with the family-mold processes, in which multiple cavities with different geometries are molded simultaneously, leading to complex flow interactions and quality variations among products. The mold used in this study consists of multiple cavities with distinct structural features and dimensions, as shown in Figure 1. These geometric differences result in non-uniform filling behavior and uneven shrinkage, making quality control particularly challenging.

Figure 1.

2D dimensional drawings of the molded part.

Based on preliminary experimental trials, feasible operating ranges of key process parameters are identified to ensure stable molding conditions and defect-free production. The investigated parameters include melt temperature, packing time, packing pressure, and cooling time, as summarized in Table 1, where the values highlighted in bold indicate the reference conditions corresponding to successful trial molding. These variables are selected as the primary drivers of the thermo-mechanical history affecting warpage and shear stress. Injection pressure is maintained constant to ensure filling stability, as preliminary trials indicated that even minor fluctuations could trigger short shots or flash, while mold temperature is fixed due to tooling constraints, specifically the lack of embedded thermal sensors within the mold cavity for precise real-time surface temperature feedback. Other process parameters, such as filling time, are also maintained at constant values throughout the experiments to isolate the effects of the selected factors. The experimentally validated parameter ranges serve as practical constraints for the subsequent numerical simulations and optimization strategy.

Table 1.

Key process parameters and corresponding ranges.

Code	Parameter	Range
$A$	Melt temperature	$230-270[\mathbf{250}]$ ( $\circ$ C)
$B$	Packing time	$4-8[\mathbf{6}]$ ( $\mathrm{s}$ )
$C$	Packing pressure	$20-40[\mathbf{30}]$ ( $\mathrm{MPa}$ )
$D$	Cooling time	$10-30[\mathbf{20}]$ ( $\mathrm{s}$ )

The material employed in this study is WONDERLOY^® PC- $345$ , a polycarbonate (PC) widely used in industrial injection molding applications owing to its high impact strength, good dimensional stability, and excellent thermal resistance. Relevant rheological, thermal, and mechanical properties of the material are incorporated into the simulation model to ensure an accurate representation of the molding behavior under various processing conditions.

Within the defined process parameter ranges and material constraints, the objective of this study is to simultaneously minimize two competing quality criteria, namely warpage and shear stress, based on numerical simulation results. Minimization of warpage is intended to improve dimensional stability and product accuracy, while reduction of shear stress aims to alleviate stress-induced defects and enhance the structural integrity of molded parts. Meanwhile, product weight is imposed as a manufacturing constraint to ensure material efficiency and cost-effectiveness. The resulting problem is therefore formulated as a simulation-based constrained multi-criteria optimization task aimed at identifying optimal process parameter combinations that achieve a balanced trade-off between dimensional accuracy and stress reduction under the prescribed weight constraint.

Experiment setup

Preliminary experiments are conducted to obtain an initial understanding of the injection molding process and to identify appropriate ranges for key process parameters. Trial samples are manufactured using a Haitian $\mathrm{PL}860/260\mathrm{j}$ injection molding machine, with a maximum injection rate of $98\mathrm{g}/\mathrm{s}$ and a clamping force of $860\mathrm{kN}$ , under different processing conditions to evaluate process stability and part quality. Figure 2 presents the experimental setup and the mold configuration used in the preliminary trials. For each trial condition, five injection shots are performed to ensure process stability before evaluation. Based on the observed molding behavior and the manufacturer’s material recommendations, feasible parameter ranges are determined and subsequently used as inputs for the numerical simulation. These parameter ranges serve as the design space for the subsequent optimization study.

Figure 2.

Experimental setup and mold details for the preliminary injection molding trials: (a) injection molding machine, (b) mold cavity, and (c) cooling channel configuration.

A computational fluid dynamics approach

All numerical simulations are performed on a personal computer running Windows $11$ ( $64$ -bit), equipped with an Intel Core i $7$ - $13700$ processor with $16$ logical processing cores available for computation, of which 8 cores are utilized for each simulation case. The system is configured with $32$ GB of RAM, providing sufficient memory capacity to handle high-fidelity numerical simulations and large-scale data processing.

Simulation setup

The filling, packing, cooling, and warpage processes in injection molding are modeled under several simplifying assumptions. The polymer melt is treated as an incompressible viscous fluid and is assumed to behave as a generalized Newtonian fluid under non-isothermal flow conditions, while inertial and gravitational effects are neglected during the filling and packing stages. The cooling process is described using a cycle-averaged approach. For the warpage analysis, the molded part is assumed to behave as a linear elastic solid undergoing small strains, and the structural response is considered quasi-static.

Filling and packing stages

The filling and packing stages are governed by the conservation laws of mass, momentum, and energy. The continuity equation is expressed as

$$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \mathbf{u})=0.$$

(1)

The momentum conservation equation

$$\frac{\mathrm{\partial }(\rho \mathbf{u})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \mathbf{u}\mathbf{u}-\mathit{\sigma })=\rho \mathbf{g},$$

(2)

$$\mathit{\sigma }=-p\mathbf{I}+\eta (\mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }{\mathbf{u}}^{\mathsf{T}}).$$

(3)

The energy conservation equation

$$\rho C_{\mathrm{p}}(\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+\mathbf{u}\cdot \mathrm{\nabla }T)=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }T\right)+\eta {\dot{\gamma }}^2,$$

(4)

where $\mathbf{u}$ denotes the velocity vector, $T$ is the temperature, $t$ represents time, $p$ is the pressure, $\mathit{\sigma }$ denotes the total stress tensor, $\rho$ is the material density, $\eta$ represents the viscosity, $k$ is the thermal conductivity, $C_{\mathrm{p}}$ denotes the specific heat capacity, and $\dot{\gamma }$ is the shear rate.

Viscosity model

The viscosity behavior of the polymer melt is described using a modified-Cross model combined with the Williams-Landel-Ferry (WLF) equation to account for temperature and pressure dependence

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

(5)

with the zero-shear viscosity defined as

$${\eta }_0(T)=D_1\exp [-\frac{A_1(T-T_{\mathrm{c}})}{A_2+(T-T_{\mathrm{c}})}],$$

(6)

$$T_{\mathrm{c}}=D_2,$$

(7)

$$A_2={\tilde{A}}_2,$$

(8)

where $n$ is the shear-thinning index, and ${\tau }^{*}$ is the critical shear stress. The parameters $D_1$ , $A_1$ , $A_2$ , and $T_{\mathrm{c}}$ are material constants in the WLF equation, with $T_{\mathrm{c}}$ and $A_2$ treated as constant material parameters, as listed in Table 2.

Table 2.

Material parameters used in the viscosity model.

Parameter	Value	Unit
$n$	$0.206834$	$-$
${\tau }^{*}$	$1.66683\times {10}^6$	$\mathrm{dyne}/\mathrm{c}{\mathrm{m}}^2$
$D_1$	$5.31058\times {10}^{18}$	$\mathrm{g}/(\mathrm{cm}\cdot \mathrm{s})$
$D_2$	$373.15$	$\mathrm{K}$
$A_1$	$45.5342$	$-$
${\tilde{A}}_2$	$51.6$	$\mathrm{K}$

Melt front tracking

A volume fraction function $f$ is used to track the melt front, where $f=0$ and $f=1$ correspond to the air and polymer melt phases, respectively, and cells with $0<f<1$ indicate the melt front. The evolution of $f$ is governed by the following transport equation

$$\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}f)=0.$$

(9)

Cooling stage

For the cycle-average mold temperature approach, the steady-state Laplace equation is used

$$k_{\mathrm{m}}(\frac{{\mathrm{\partial }}^2\overline{T}}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2\overline{T}}{\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2\overline{T}}{\mathrm{\partial }{z}^2})=0,$$

(10)

where $\overline{T}$ denotes the cycle-averaged mold temperature, $k_{\mathrm{m}}$ is the thermal conductivity of the mold base, and $x$ , $y$ , and $z$ represent the Cartesian coordinates.

Warpage analysis

Warpage is analyzed using solid mechanics theory. The governing equilibrium equation describing the warpage behavior is expressed as

$${\sigma }_{\mathrm{ij},\mathrm{j}}+f_{\mathrm{i}}=0,$$

(11)

where ${\sigma }_{\mathrm{ij}}$ denotes the stress tensor and $f_{\mathrm{i}}$ represents the body force.

Mesh generation strategy

Figure 3 illustrates the computational mesh employed in the CFD simulation of the injection molding process. The overall mesh configuration of the mold cavity and runner system is shown in Figure 3(a), while Figure 3(b) highlights the locally refined mesh applied in regions with complex geometry and expected high gradients, such as the junction between the runner and gate. A detailed view of the mesh transition region between different refinement levels is presented in Figure 3(c).

Figure 3.

Mesh configuration used in the CFD simulation: (a) global mesh, (b) locally refined mesh, and (c) mesh transition region.

The mesh is designed to achieve a balance between computational cost and numerical accuracy. Local mesh refinement is applied near the runner, gate, and cavity to better capture the evolution of the melt front and local flow features, while a coarser mesh is used in regions away from these critical zones. Smooth mesh transitions are ensured to avoid abrupt changes in element size, which helps maintain numerical stability and reduce discretization errors. This mesh strategy is considered adequate for accurately resolving the flow and thermal fields during the injection molding simulation.

Surrogate-based optimization framework

An integrated surrogate-based optimization framework is proposed to address the multi-criteria optimization problem in family-mold processes. The overall workflow of the proposed approach is illustrated in Figure 4 as a flowchart, highlighting the main stages from simulation-based data generation and surrogate modeling to optimization and solution validation.

Figure 4.

Flowchart of the optimization pipeline.

Face-centered central composite design

Design of Experiments (DOE) provides a systematic and efficient framework for generating representative datasets from a multidimensional design space while minimizing the number of experimental trials. Among DOE strategies, central composite design (CCD) is a standard and efficient approach that augments a factorial or fractional factorial design with additional design points, enabling interaction and curvature effects to be captured using a limited number of experiments. ²¹ The specific CCD variant is defined by the placement of the axial (star) points relative to the factorial region. ²² In this work, a face-centered central composite design (FCCCD) is employed, where the axial points coincide with the centers of the factorial faces, leading to three discrete levels for each design variable. This configuration provides a well-distributed and bounded sampling of the design space without extending beyond the predefined experimental limits. The FCCCD-generated dataset is subsequently used as the DOE foundation for surrogate modeling via Kriging and for multi-criteria optimization using the NSGA-II algorithm.

Kriging-based surrogate modeling

Kriging, originally developed in the field of geostatistics, is an interpolation-based surrogate modeling technique that has been widely adopted in engineering design and optimization. ²³ It assumes the unknown response function to be a realization of a stochastic process, thereby providing both a mean prediction and an explicit quantification of prediction uncertainty. From a probabilistic perspective, Kriging can be interpreted as Gaussian process regression, in which the response at each design point is modeled as a Gaussian random variable, and the correlation between design points is characterized by a covariance (kernel) function. The choice of kernel function and its hyperparameters governs the smoothness and nonlinear representational capability of the surrogate, making Kriging a data-efficient and widely used approach for approximating computationally expensive high-fidelity models in surrogate-assisted and uncertainty-aware optimization problems.^24,25 For example, Yang et al. ²⁶ proposed a Kriging-assisted NSGA-II framework for multi-objective optimization of laser–magnetic hybrid welding parameters, demonstrating effective defect reduction and improved welding quality with high computational efficiency. Similarly, Baum et al. ²⁷ compared Kriging and response surface methodology surrogate models for injection molding process optimization and found that Kriging provides superior prediction accuracy and efficiency, especially in complex geometrical configurations. Additionally, Gao and Wang ²⁸ developed an adaptive Kriging-based surrogate optimization method with an expected improvement strategy for injection molding, achieving effective warpage reduction with substantially reduced computational effort. Collectively, these findings motivate the adoption of a Kriging-based surrogate modeling strategy in the present work to efficiently optimize the process parameters under the considered system configuration.

Leave-one-out cross-validation

Cross-validation is commonly employed to evaluate the predictive performance and generalization capability of surrogate models under data partitioning schemes.^29,30 Among various cross-validation techniques, leave-one-out cross-validation (LOO-CV) is adopted in this study due to its rigorous and unbiased estimation of prediction accuracy, particularly when the available dataset is limited in size. LOO-CV can be regarded as a special case of the $k$ -fold cross-validation scheme, in which the number of folds is equal to the number of available samples ( $k=p$ ). ³¹ In this procedure, one data point is sequentially excluded from the dataset and used as a validation sample, while the remaining data are used for model construction, and this process is repeated until each data point has been used once for validation. Unlike the coefficient of determination obtained from model fitting $R_{\mathrm{fit}}^2$ , which reflects only the goodness of fit to the training data, the cross-validated coefficient of determination $R_{\mathrm{cv}}^2$ derived from LOO-CV provides a more reliable measure of the true predictive capability of the model. In addition, the root mean square error $RMSE$ is used to quantify the overall prediction error, serving as a complementary indicator of model accuracy. ³²

Global sensitivity analysis

To quantify the individual and coupled influences of the design parameters on the molding quality criteria, a robust mathematical tool is required to scan the multi-dimensional design space. In this study, the variance-based Sobol ³³ sensitivity analysis is deployed to decompose the total response variance into fractional contributions from single variables and their higher-order combinations. The significance of each process parameter is formally evaluated through two distinct metrics: the first-order sensitivity index ( $S_i$ ) and the total-effect sensitivity index ( $S_{Ti}$ ). While $S_i$ measures the main effect of an isolated parameter on the response variance without considering cross-variable influences, $S_{Ti}$ accounts for the comprehensive contribution of the parameter, including all higher-order interactive effects. ³⁴ By computing both $S_i$ and $S_{Ti}$ across the defined parameter ranges, the framework systematically establishes a rigorous mathematical mapping of factor significance prior to the multi-criteria optimization phase.

Non-dominated Sorting Genetic Algorithm II

The Non-dominated Sorting Genetic Algorithm II (NSGA-II) ³⁵ is an advanced evolutionary algorithm developed for solving multi-objective optimization problems. It employs fast non-dominated sorting, elitism, and crowding distance mechanisms to efficiently obtain a diverse set of Pareto-optimal solutions. Due to its robustness and computational efficiency, NSGA-II has been widely applied to complex engineering optimization problems. For instance, Banh et al. ¹² employed NSGA-II to optimize multiple conflicting process parameters in a complex engineering system, aiming to simultaneously improve several performance objectives. By exploring the trade-offs among competing objectives, the optimization results yielded a well-distributed Pareto-optimal solution set, suggesting potential improvements in overall system performance. Alvarado-Iniesta et al. ³⁶ integrated an artificial neural network with NSGA-II to optimize injection molding process parameters under multiple objectives related to quality, productivity, and energy efficiency. Likewise, Tian et al. ³⁷ applied RSM-based regression models combined with NSGA-II to identify suitable injection molding process parameters under multiple quality and energy-related objectives. The results indicate that the obtained parameter combinations have the potential to improve process stability, contribute to better compliance of product length with specifications, and reduce product weight and energy consumption, supporting greener manufacturing practices. From these studies, NSGA-II is demonstrated as an effective multi-objective optimization tool for handling conflicting process parameters, enabling the simultaneous consideration of multiple performance, quality, and energy-related objectives in complex engineering and manufacturing systems. Accordingly, in this study, NSGA-II is employed for multi-criteria optimization with the population size and number of generations set to $100$ and $150$ , respectively, to ensure sufficient exploration of the solution space and convergence toward optimal solutions.

Results and discussion

Preliminary experimental results

Preliminary experimental trials are conducted to evaluate the practical feasibility of the family-mold processes and to assess mold performance under real manufacturing conditions. During this phase, the mold is installed on the injection molding machine, and multiple trial shots are carried out to investigate the filling behavior, part formation, and process stability. Figure 5 presents representative photographs captured during the experimental trials, illustrating typical process-induced defects identified in the preliminary investigation, such as short shot, flash, and overpacking phenomena.

Figure 5.

Defects in family mold process from left to right: short shot, flash, overpacking.

Based on the preliminary experimental results, feasible ranges of key process parameters, including melt temperature, packing time, packing pressure, and cooling time, are identified and used to define the input bounds for numerical simulations and multi-criteria optimization. To further validate the applicability of these parameter ranges, a successful trial molding is conducted using a selected parameter set, with the corresponding values highlighted in bold in Table 1. The molded parts obtained from a single representative family-mold injection cycle under these conditions are shown in Figure 6, demonstrating stable mold operation and acceptable part quality without apparent process-induced defects. The total weight of the molded parts produced in this injection cycle is measured as $19.220\mathrm{g}$ using a precision electronic balance. By restricting the optimization process to experimentally achievable operating conditions, practical relevance is ensured while computational effort and experimental cost are significantly reduced. Consequently, the proposed approach enables efficient exploration of the design space and provides a robust foundation for subsequent data-driven optimization.

Figure 6.

The molded parts meet quality requirements from left: molded parts and corresponding weight.

Mesh independence study

A mesh independence study is performed to examine the effect of mesh density on the numerical accuracy of the simulation results and to confirm that the solution is independent of the spatial discretization. To this end, five different mesh densities are investigated, consisting of $263,000$ , $485,000$ , $795,000$ , $1,080,000$ , and $1,327,000$ nodes. For each mesh configuration, the same processing conditions and simulation settings are applied.

As shown in Figure 7, the difference between the simulation results obtained with $1,080,000$ and $1,327,000$ nodes is minimal, indicating that the numerical solution has essentially converged at this level of mesh refinement. The variation in warpage is less than $0.003\mathrm{mm}$ , while the change in volumetric shrinkage remains below $0.85\mathrm{\%}$ . In addition, the computed weight shows only a marginal difference of less than $0.01\mathrm{g}$ between these two mesh densities. In contrast, the computational time increases significantly from approximately $7.0$ hours to $9.2$ hours when the mesh is further refined. Considering that the numerical error between these two mesh densities is negligible (error $<1\mathrm{\%}$ for all monitored responses), whereas the computational cost increases substantially, a mesh with $1,080,000$ nodes is selected as the optimal mesh for all subsequent simulations.

Figure 7.

Mesh independence study results.

Simulation-assisted data acquisition

Following the identification of feasible ranges for key process parameters through preliminary experimental trials, a simulation-based data acquisition strategy is employed to construct a comprehensive dataset for surrogate modeling and optimization, in which FCCCD is adopted to systematically sample the multi-dimensional design space within the experimentally validated parameter bounds, and the resulting DOE matrix is summarized in Table 3. High-fidelity numerical simulations are then performed using Moldex3D software for each design point to predict the corresponding quality responses, including warpage and shear stress, such that the resulting simulation outputs, together with their associated process parameter settings, form a structured dataset capturing the nonlinear and coupled relationships between process variables and molding quality and thereby providing a reliable foundation for subsequent stages of the proposed methodology.

Table 3.

Experimental design matrix and corresponding responses.

$Run$	$A$ ( $\circ$ C)	$B$ ( $\mathrm{s}$ )	$C$ ( $\mathrm{MPa}$ )	$D$ ( $\mathrm{s}$ )	$\mathrm{Warpage}$ ( $\mathrm{mm}$ )	$\mathrm{Shear}\mathrm{stress}$ ( $\mathrm{MPa}$ )	$\mathrm{Weight}$ ( $\mathrm{g}$ )
1	250	6	30	20	0.312	0.387	19.203
2	230	4	40	30	0.296	0.118	19.276
3	250	4	30	20	0.343	0.146	19.178
4	270	6	30	20	0.302	0.208	19.246
5	250	6	20	20	0.343	0.550	19.116
6	250	6	30	20	0.312	0.387	19.203
7	230	4	40	10	0.299	0.108	19.262
8	250	6	30	10	0.314	0.310	19.200
9	250	6	30	20	0.312	0.387	19.203
10	250	6	30	20	0.312	0.387	19.203
11	270	8	20	10	0.305	0.608	19.166
12	270	4	20	10	0.385	0.112	19.084
13	230	6	30	20	0.307	0.418	19.181
14	250	6	40	20	0.267	0.200	19.306
15	250	6	30	20	0.312	0.387	19.203
16	250	6	30	20	0.312	0.387	19.203
17	270	8	40	30	0.221	0.279	19.362
18	230	8	20	10	0.289	0.883	19.156
19	250	6	30	30	0.310	0.429	19.207
20	230	8	40	30	0.227	0.495	19.288
21	270	4	40	10	0.300	0.107	19.290
22	270	4	20	30	0.384	0.138	19.104
23	230	4	20	10	0.387	0.233	19.053
24	270	8	40	10	0.230	0.216	19.354
25	230	4	20	30	0.381	0.265	19.076
26	230	8	20	30	0.281	0.962	19.162
27	230	8	40	10	0.240	0.418	19.276
28	270	8	20	30	0.290	0.736	19.184
29	270	4	40	30	0.298	0.130	19.310
30	250	8	30	20	0.271	0.737	19.210

Figure 8 presents representative distributions of warpage and shear stress obtained from Moldex3D simulations. Warpage is found to be most pronounced in the vicinity of the gate region and near the end of the flow path, while comparatively lower deformation is observed in the mid-length regions of the molded parts. This behavior can be attributed to the combined influence of non-uniform pressure decay, temperature gradients, and differential shrinkage between the early-filled and late-filled regions during the packing and cooling stages. In contrast, shear stress exhibits elevated levels along the flow direction, particularly near the gate and in regions experiencing high melt velocity and strong velocity gradients, whereas relatively lower shear stress is observed in regions farther from the gate where flow deceleration and thermal relaxation occur. These spatial characteristics provide qualitative validation of the physical consistency of the numerical model and reveal the coupled evolution of stress accumulation and deformation in family product molding. Moreover, the cavity-to-cavity variations highlight the influence of flow imbalance and thermal interactions, further emphasizing the necessity of a data-driven surrogate modeling approach to accurately capture the nonlinear relationships between process parameters and quality responses.

Figure 8.

Representative Moldex3D simulation results: (a) warpage, and (b) shear stress.

Surrogate model performance

A qualitative assessment of the surrogate model performance is first conducted using predicted-versus-actual plots. Figure 9 presents LOO-CV results for warpage and shear stress, where the predicted values are plotted against the corresponding experimental observations. The data points are closely distributed along the ideal diagonal line, indicating a strong agreement between the model predictions and actual values and qualitatively demonstrating the robustness and reliability of the developed Kriging surrogate models.

Figure 9.

Predicted $-$ actual comparison of warpage and shear stress using Kriging surrogate models.

In addition to the qualitative evaluation, a quantitative assessment of the predictive performance of the Kriging surrogate models is performed using LOO-CV for both warpage and shear stress responses. Table 4 summarizes the resulting $RMSE$ and $R_{\mathrm{cv}}^2$ values. The $RMSE$ values are small relative to the corresponding ranges of the response variables, while the $R_{\mathrm{cv}}^2$ values remain consistently close to unity, indicating that the Kriging models can accurately capture and explain the variance in both warpage and shear stress responses. These results confirm that the developed surrogate models provide sufficient predictive accuracy to support the subsequent multi-criteria optimization.

Table 4.

Prediction accuracy of the kriging surrogate model.

Response	$RMSE$	$R_{\mathrm{cv}}^2$
Warpage	0.007485	0.9696
Shear stress	0.033147	0.9787

Parameter sensitivity analysis

As presented in Figure 10(a) and (b), packing time ( $B$ ) and packing pressure ( $C$ ) have emerged as the primary factors influencing both criteria, whereas melt temperature ( $A$ ) and cooling time ( $D$ ) have exhibited marginal primary and interactive effects. Specifically, $B$ has exerted the most pronounced impact on shear stress with an index approaching $0.7$ , followed by its considerable governance over warpage ( $\approx 0.51$ ). This phenomenon is highly consistent with the rheological solidification behavior of the polymer melt during the post-filling stage. Prolonged packing phase directly extends the duration under which the melt is subjected to forced flow through the restricted gates, thereby promoting the development of frozen-in molecular orientation and residual stresses near the solidified skin layer. Concurrently, $C$ has demonstrated a distinct sensitivity profile, contributing substantially to warpage ( $\approx 0.48$ ) but showing a more moderate effect on shear stress ( $\approx 0.21$ ). From a thermodynamic perspective, $C$ dictates the localized volumetric shrinkage by forcing additional material to compensate for thermal contraction. An inadequate packing pressure fails to maintain a uniform specific-volume distribution across the tubular domains, leading to severe differential shrinkage and, consequently, elevated warpage.

Figure 10.

Global sensitivity analysis via Sobol indices: (a) first-order indices ( $S_i$ ), and (b) total-effect indices ( $S_{Ti}$ ).

Crucially, a comparative analysis between $S_i$ and $S_{Ti}$ reveals that the total-effect indices for all parameters are nearly identical to their first-order counterparts. This quantitative observation indicates that the multi-variable interaction effects within the defined operational bounds are relatively minimal. In physical terms, the close agreement between $S_i$ and $S_{Ti}$ suggests that the dominant thermodynamic effect of volumetric compensation (primarily governed by $C$ ) and the rheological effect associated with shear-induced residual stresses (primarily governed by $B$ ) remain only weakly coupled under the investigated processing conditions. These numerical insights, enabled by the efficient Kriging responses, indicate a predominantly additive response structure with limited higher-order interactions, thereby providing a mathematically consistent basis for the subsequent multi-criteria design optimization phase.

Multi-criteria optimization results

Based on the reliable Kriging surrogate models developed for the considered quality responses, the constrained multi-criteria optimization problem is formulated as

$$\{\begin{array}{c}min{f}_1(\mathbf{x})=\text{Warpage}(\mathbf{x}),\\ min{f}_2(\mathbf{x})=\text{Shear,,stress}(\mathbf{x}),\\ g_1(\mathbf{x})=\text{Weight}(\mathbf{x})-{\text{Weight}}_{max}\le 0.\end{array}$$

(12)

where $\mathbf{x}$ denotes the vector of injection molding process parameters, $f_1$ and $f_2$ represent the objective functions corresponding to warpage and shear stress, respectively, and $g_1$ defines the inequality constraint imposed on the maximum allowable part weight, ${\text{Weight}}_{max}=19.180\mathrm{g}$ . The formulated multi-objective optimization problem is solved using the NSGA-II algorithm to systematically explore the trade-off between competing quality criteria in the family-mold process. Increasing the packing pressure introduces additional material into the cavity to compensate for thermal shrinkage, thereby reducing warpage. However, excessive packing during the cooling stage increases shear rates in the polymer melt, leading to higher molecular orientation that becomes frozen upon solidification and consequently elevates shear stress. This antagonistic behavior establishes an inherent trade-off between warpage reduction and shear stress amplification. The resulting Pareto-optimal solutions, presented in Figure 11, delineate the feasible design space and clearly highlight the conflicting nature of the considered criteria.

Figure 11.

Pareto-optimal solutions showing the trade-off between warpage and shear stress: (a) full Pareto front and (b) selected Pareto solutions.

As shown in Figure 11(a), a distinct and continuous Pareto front is obtained, indicating a clear and pronounced trade-off between warpage and shear stress. In general, lower warpage levels are associated with higher shear stress, reflecting the inherent conflict between dimensional accuracy and material efficiency. To facilitate practical decision-making, a subset of representative Pareto-optimal solutions is highlighted in Figure 11(b), which serves as candidate process parameter combinations for subsequent verification and evaluation.

To investigate the reliability of the Pareto-optimal solutions obtained from the surrogate-assisted multi-criteria optimization strategy, selected candidate solutions are further evaluated using high-fidelity CFD simulations, and the comparison results are summarized in Table 5. The verification results demonstrate a close agreement between the surrogate-based predictions and the CFD simulation outcomes for both warpage and shear stress. For warpage, the minimum and maximum relative errors are $1.00\mathrm{\%}$ and $2.07\mathrm{\%}$ , respectively, with an average relative error of $1.56\mathrm{\%}$ . For shear stress, the corresponding minimum and maximum relative errors are $0.61\mathrm{\%}$ and $8.44\mathrm{\%}$ , respectively, and the average relative error is $3.58\mathrm{\%}$ . Despite the inherent sensitivity of warpage to variations in flow behavior, packing pressure, and cooling conditions, and the observable variability in shear stress induced by changes in local flow and stress distributions, all discrepancies remain well controlled. These results confirm that the proposed surrogate-assisted optimization framework can reliably identify high-quality Pareto-optimal solutions and demonstrate the effectiveness and robustness of NSGA-II in addressing multi-criteria optimization problems in family-mold processes.

Table 5.

Prediction and simulation results with corresponding errors under the full case weight constraint of $19.180\mathrm{g}$ .

$\mathbf{C}\mathbf{a}\mathbf{s}\mathbf{e}$	$A$	$B$	$C$	$D$	$\mathrm{Warpag}{\mathrm{e}}_{\mathrm{pred}}$	$\mathrm{Warpag}{\mathrm{e}}_{\mathrm{sim}}$	$\mathrm{RE}$	$\mathrm{Shear}\mathrm{stres}{\mathrm{s}}_{\mathrm{pred}}$	$\mathrm{Shear}\mathrm{stres}{\mathrm{s}}_{\mathrm{sim}}$	$\mathrm{RE}$
	( $\circ$ C)	( $\mathrm{s}$ )	( $\mathrm{MPa}$ )	( $\mathrm{s}$ )	( $\mathrm{mm}$ )	( $\mathrm{mm}$ )	( $\mathrm{\%}$ )	( $\mathrm{MPa}$ )	( $\mathrm{MPa}$ )	( $\mathrm{\%}$ )
1	230.00	6.52	29.74	10.04	0.299	0.302	1.00	0.449	0.454	1.11
2	230.00	6.67	29.66	11.22	0.295	0.300	1.69	0.489	0.492	0.61
3	230.05	6.58	29.74	10.00	0.297	0.302	1.68	0.461	0.491	6.51
4	230.00	6.96	29.51	10.02	0.290	0.296	2.07	0.533	0.488	8.44
5	230.00	6.77	29.64	10.01	0.294	0.298	1.36	0.495	0.501	1.21
$\mathbf{A}\mathbf{v}\mathbf{e}\mathbf{r}\mathbf{a}\mathbf{g}\mathbf{e}$	$-$	$-$	$-$	$-$	$-$	$-$	1.56	$-$	$-$	3.58

$*$ pred: prediction, sim: simulation, and RE: relative error.

Limitations and outlook

The current study primarily addresses a multi-criteria optimization problem involving two competing quality criteria, warpage and shear stress, subject to a product weight constraint. However, industrial injection molding practice is inherently a many-criteria optimization task. Future enhancements should extend this framework to incorporate a broader range of critical performance metrics, such as cycle time, sink marks, volumetric shrinkage, and energy efficiency. Integrating these additional objectives will provide a more holistic and industrially relevant production strategy for complex manufacturing environments.

Another limitation pertains to the material specificity and the geometric layout of the mold used for validation. While the proposed digital twin architecture demonstrates strong predictive capabilities for amorphous Polycarbonate (PC), its generalizability to other material classes, particularly semi-crystalline polymers, requires further validation. Unlike amorphous polymers, semi-crystalline materials exhibit complex crystallization kinetics and latent heat release during the cooling phase, which could introduce additional non-linearities into the optimization landscape. Specifically, the microstructural evolution involves highly non-linear phase transitions where local shear layers alter the thermodynamic conditions, triggering simultaneous nucleation and crystal growth that release substantial latent heat and cause abrupt volumetric shrinkage. To adapt the proposed framework to these materials, these rheological and thermal non-linearities can be completely handled at the modular CFD layer by utilizing proper multi-phase state equations to capture the localized crystallinity and thermal disruptions. Because the Kriging surrogate acts as a robust non-linear mathematical interpolator, it is fully capable of mapping the resulting complex output responses passed from the CFD layer without requiring any structural modifications to the core optimization pipeline. Furthermore, while the framework successfully captured the coupled interactions in a family mold with relatively straightforward geometries, real-world applications often involve high-precision industrial components with intricate features. Future research will focus on expanding the framework’s modularity to accommodate different rheological and thermodynamic material behaviors, as well as handling more complex multi-cavity layouts that induce severe flow and thermal imbalances. Such advancements are necessary to further demonstrate the robustness and adaptability of the proposed optimization framework across diverse manufacturing scenarios.

Additionally, the current framework relies on a deterministic approach, assuming constant process inputs. Because industrial molding is subject to stochastic variations (e.g., fluctuations in melt temperature or packing pressure), future iterations will utilize the variance-prediction capabilities of the Kriging model to execute formal stochastic uncertainty quantification, ensuring the robustness of the Pareto-optimal solutions under real-world process noise.

Finally, the reliance on a static FCCCD for data acquisition poses efficiency constraints. As a static approach, it may not capture localized nonlinearities across the design space as effectively as dynamic methods. Prospective iterations of the framework should implement adaptive sequential sampling strategies to selectively enrich high-uncertainty regions of the design space. This would significantly improve surrogate accuracy while further reducing the total number of high-fidelity CFD simulations required, thereby enhancing the overall computational efficiency.

Conclusion

In this study, a machine learning-driven constrained multi-criteria optimization framework is proposed for the intelligent design of family-mold processes, in which quality control is complicated by strong and complex cavity interactions. The framework integrates a face-centered central composite design (FCCCD), a Kriging-based surrogate model, and the Non-dominated Sorting Genetic Algorithm II (NSGA-II) to simultaneously minimize two competing quality criteria, namely warpage and shear stress, under a prescribed product weight constraint. High-fidelity simulation data generated using Moldex3D computational fluid dynamics (CFD) analysis are employed to train the surrogate model, enabling accurate representation of the nonlinear relationships between process parameters and quality responses. The primary innovation of this work lies in its formulation of a constrained multi-criteria approach specifically tailored for family-mold systems. While existing studies often simplify these challenges into single-objective tasks or unconstrained scenarios, the proposed methodology explicitly addresses the competitive trade-offs between dissimilar cavities. By simultaneously optimizing multiple conflicting quality metrics under a strict product weight constraint, the framework successfully resolves complex inter-cavity flow imbalances while ensuring material efficiency and practical viability. The resulting Pareto-optimal solutions provide an effective trade-off between dimensional stability and stress reduction while ensuring material efficiency, and validation results further confirm the high predictive accuracy and robustness of the proposed surrogate modeling and optimization framework.

In conclusion, the proposed data-driven optimization framework offers an efficient and reliable approach for quality-oriented optimization of family-mold processes. By significantly reducing computational cost while maintaining optimization accuracy, the framework supports informed process decision-making and enhances the feasibility of implementing intelligent optimization strategies in industrial practice. The proposed methodology therefore represents a practical contribution toward smart manufacturing systems and can be readily extended to other complex many-objective optimization problems in advanced manufacturing processes.