Bruno A Storti
email: bruno.storti@univ-nantes.fr

Vincent Sobotka

A numerical framework for three-dimensional optimization of cooling channels in thermoplastic printed molds

Keywords: Cooling channels, Overset finite elements, Injection molding process, Thermoplastic mold, Transient heat transfer, Stochastic optimization

The rapid progress of additive manufacturing (AM) technology in the last few decades has unlocked new potentials for the injection molding industry by enabling the quick and cost-efficient manufacture of molds with complex geometries. Despite the wide literature regarding experimental investigations on AM polymer-based soft tools specifically produced to cope with low-volume production of small-sized injected parts, their lifespan remains uncertain, and premature failures are often related to poor thermal performance. This paper is devoted to the numerical study and enhancement of the heat transfer within a thermoplastic 3D printed insert with cooling channels (CCs). Experimental thermal characterization is performed on printed composite samples made of polycarbonate reinforced with carbon fibers, considered as tool material. The simulations are performed on an industrial case, thereby facilitating a comprehensive validation of the proposed framework. Parametric studies show a marked cycle time sensitivity to insert thermal conductivity within the 0.1-1.0W/(mK) range while demonstrating negligible influence on cycle times for polymer-polymer thermal contact resistance values below 10 -3 m 2 K/(W). Furthermore, to find a suitable arrangement of the CCs' layout, we propose here an accurate optimization methodology based on 3D overset meshes in the finite element method context coupled to the augmented Lagrangian particle swarm optimizer. The optimized bent CCs' configuration enhances part surface temperature uniformity by 42% and reduces its temperature delta by over 6 • C, all while employing 67% of the reference cycle time. The thermal shortcomings of the thermoplastic AM mold, compared to its steel counterpart, are also addressed in this work.

Introduction

The injection molding process, comprising the injection of the molten polymer into a mold's cavity to solidify and yield the final product, is a prevalent method in the thermoplastic manufacturing industry for its low cost and high production efficiency. The production of steel molds typically requires multiple machining operations, such as drilling and milling, which are time-consuming and costly. Consequently, their use is generally restricted to applications that involve high-volume manufacturing. Due to the change of standard on the actual manufacturing industry from mass production to a made-to-order perspective, and the necessity of the enterprises for quick development of new products, injection mold (IM) manufacturers are often demanded for rapid tooling to deal with low-volume batches and rapid prototyping. In this context, additive manufacturing (AM) emerged as a promising technology that enables to produce polymer-based molds with complex geometries at a lower cost while reducing the time of manufacture. Notwithstanding the advantages of utilizing AM polymer-based molds, their thermomechanical properties are inferior to those of steel molds. Specifically, they exhibit thermal diffusivity and elasticity modulus on the order of 100 times lower, and a thermal expansion coefficient approximately 10 times higher, making the cooling of the part in a reasonable time without mold failure a challenge. Indeed, the unpredictability of the tool life of such polymer-based inserts is the main concern of mold manufacturers nowadays, and the reason of many ongoing research [START_REF] Davoudinejad | Experimental investigation and thermomechanical modelling for tool life evaluation of photopolymer additively manufactured mould inserts in different injection moulding conditions[END_REF][START_REF] Whatcott | Effects of Conformal Cooling Channels on Additively Manufactured Injection Molding Tooling[END_REF][START_REF] Bogaerts | Influence of thermo-mechanical loads on the lifetime of plastic inserts for injection moulds produced via additive manufacturing[END_REF][START_REF] Lozano | Analysis and advances in additive manufacturing as a new technology to make polymer injection molds for world-class production systems[END_REF].

Rapid tooling, a term denoting the expedited production of molds or tooling components while bypassing intricate and time-consuming mechanical procedures, can be classified into either 'soft' or 'hard' categories, contingent upon the base material utilized, such as high-strength steel [START_REF] Mendible | Comparative study of rapid and conventional tooling for plastics injection molding[END_REF][START_REF] Kanbur | Metal additive manufacturing of plastic injection molds with conformal cooling channels[END_REF] or softer materials like polymers [START_REF] Davoudinejad | Experimental investigation and thermomechanical modelling for tool life evaluation of photopolymer additively manufactured mould inserts in different injection moulding conditions[END_REF][START_REF] Whatcott | Effects of Conformal Cooling Channels on Additively Manufactured Injection Molding Tooling[END_REF][START_REF] Kuo | Characterizations of polymer injection molding tools with conformal cooling channels fabricated by direct and indirect rapid tooling technologies[END_REF], silicone rubbers [START_REF] Kuo | Improving cooling performance of injection molding tool with conformal cooling channel by adding hybrid fillers[END_REF], or epoxy resins [START_REF] Kuo | Characterizations of polymer injection molding tools with conformal cooling channels fabricated by direct and indirect rapid tooling technologies[END_REF][START_REF] Feng | Design and fabrication of conformal cooling channels in molds: Review and progress updates[END_REF][START_REF] Kuo | Development and application of a conformal cooling channel with easy removal and smooth surfaces[END_REF]. This classification also extends to the differentiation between 'direct' and 'indirect' production methods [START_REF] Kuo | Characterizations of polymer injection molding tools with conformal cooling channels fabricated by direct and indirect rapid tooling technologies[END_REF][START_REF] Rahmati | Direct Rapid Tooling[END_REF], depending on whether the mold is fabricated directly via AM, for instance, or indirectly through casting on a master pattern. Soft tooling that incorporates conformal cooling channels (CCs) can employ a hybrid manufacturing technique as well, involving the initial 3D printing of CCs using AM technology, followed by their subsequent deployment in an indirect process to cast a mold, with silicone rubber [START_REF] Kuo | Improving cooling performance of injection molding tool with conformal cooling channel by adding hybrid fillers[END_REF] or epoxy resin [START_REF] Kuo | Effects of different mold materials and coolant media on the cooling performance of epoxy-based injection molds[END_REF] as potential casting materials. While indirect tooling usually requires more manufacturing steps and specific procedures for CCs' material removal [START_REF] Kuo | Development and application of a conformal cooling channel with easy removal and smooth surfaces[END_REF], AM can overcome this by manufacturing a mold in one single step, which can eventually offer superior accuracy and thermomecanical properties [START_REF] Kuo | Improving cooling performance of injection molding tool with conformal cooling channel by adding hybrid fillers[END_REF][START_REF] Rahmati | Direct Rapid Tooling[END_REF]. Leveraging recent advancements in AM technology, particularly the utilization of charged high-performance thermoplastic polymers [START_REF] Chen | Thermal conductivity of polymer-based composites: Fundamentals and applications[END_REF][START_REF] Boparai | Development of rapid tooling using fused deposition modeling: a review[END_REF], the study and enhancement of the thermal performance of a soft polymeric mold directly manufactured by AM are the main concerns of this work.

The feasibility of using AM polymer-based molds has been explored experimentally since a few decades ago when Rahmati and Dickens [START_REF] Rahmati | Stereolithography for injection mould tooling[END_REF][START_REF] Rahmati | Rapid tooling analysis of Stereolithography injection mould tooling[END_REF] employed a stereolithography (SLA) technique to produce a resin mold that was able to produce a small batch of up to 500 parts. The produced parts were small (≈10×10×10mm 3 ), thus no CCs were used, however, the cycle time had to be significantly increased to avoid warpage and to achieve a proper mold temperature. Davoudinejad et al. [START_REF] Davoudinejad | Experimental investigation and thermomechanical modelling for tool life evaluation of photopolymer additively manufactured mould inserts in different injection moulding conditions[END_REF] performed experimental and numerical investigations of an AM mold insert made of a methacrylic photopolymer by the digital light processing (DLP) technology to produce small parts (≈20×20×2.5mm 3 ). By applying a design of experiments (DOE) technique for different combinations of process parameters, they concluded that the best repeatability was obtained with the lower mold temperature and the longest cooling times. Mendible et al. [START_REF] Mendible | Comparative study of rapid and conventional tooling for plastics injection molding[END_REF] conducted a comprehensive investigation that compared the thermomechanical behavior of rapid and conventional tooling used for plastics injection molding, both experimentally and numerically. The inserts were indirectly cooled by conformal CCs that were in contact with the steel mold plates. The polymeric insert (PI) demonstrated a limited lifespan, slightly exceeding 100 produced parts before failure occurred.

Due to the low thermal conductivity of the PI, longer cooling, holding, and ejection times were necessary to achieve the desired ejection temperature. Despite the reduced dimensions of the part, the results of the study demonstrated that a cycle time over four times longer was needed compared to a conventional steel insert. In addition, the PI demonstrated the highest degree of shrinkage variability and the most significant temperature differences. A recent study of Bagalkot et al. [START_REF] Bagalkot | Categorization of failures in polymer rapid tools used for injection molding[END_REF] categorized early failures on AM soft injection molds. They identified the high tool temperature as the main cause of the short lifetime of such molds.

Bogaerts et al. [START_REF] Bogaerts | Influence of thermo-mechanical loads on the lifetime of plastic inserts for injection moulds produced via additive manufacturing[END_REF] conducted a similar comparative study and found that temperature concentrations and large variations were the primary cause of insert failures. Additionally, numerical models were shown to be useful in designing plastic inserts and optimizing process parameters to improve thermal performance.

From such works, it becomes evident that effective thermal control measures for the design of AM polymeric molds are crucial, and accurate numerical simulations of the process represent a major gap in the literature [START_REF] Kovács | Thermal simulations and measurements for rapid tool inserts in injection molding applications[END_REF][START_REF] Schuh | Towards Temperature Control Measures for Polymer Additive Injection Molds[END_REF]. Despite the number of recent studies on the subject, some inquiries remain unresolved, for instance:

• To what extent is the cycle time impacted by the thermal conductivity of the polymer?

• Does the thermal contact resistance at the polymer-polymer interface play a crucial role during the cooling phase?

• How do CCs perform in AM polymeric inserts as compared to steel inserts?

• What level of improvement in temperature homogeneity of the part can be attained through the CCs' optimization?.

One goal of this work is to shed light on these questions. The present study involves an actual industrial scenario where a conventional steel insert with CCs is currently utilized to manufacture parts made of acrylonitrile butadiene styrene (ABS). The objective of this study is to evaluate the potential of replacing the existing steel insert with a polymer-based counterpart. The selected polymer material for this study is a polycarbonate matrix reinforced with carbon fibers. Its thermal properties are determined through experimental measurements conducted on 3D printed samples produced by collaborators of the SAMFAST project, which focuses on the development of inserts using fused deposition modeling (FDM) technology. Several 1D-3D parametric studies are performed by the finite element method (FEM) to assess the thermal performance of the CCs inside the plastic insert. Furthermore, to address the aforementioned inquiries, a numerical optimization methodology is presented, exploiting the flexibility of the AM process.

Great research efforts have already been put into developing numerical methodologies to find the best disposition of CCs inside conventional steel molds aiming to reduce warpage while improving productivity.

For instance, Agazzi et al. [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF] developed an innovative methodology using the commercial FEM code COMSOL based on a morphological approach in order to find the best temperature distribution over a dilated 3D cooling surface aiming to improve the temperature homogeneity at the part surface during the cooling stage. The methodology's main advantage is that no a priori disposition nor shape of the CCs is required. In addition, the scheme can be extended to other processes, such as hot stamping and over-molding [START_REF] Reyes | Inverse heat transfer optimization of stamping with overmolding process involving high performance thermoplastic composites: experimental validation[END_REF]. On the other hand, as the scheme relies on the gradient-based conjugated gradient algorithm, the final solution is conditioned by the initial state employed and the possibility of falling into a local minima arises. Furthermore, the CCs must be defined by the designer as a postprocess. Using the stochastic genetic algorithm (GA), Mercado-Colmenero et al. [START_REF] Mercado-Colmenero | A new conformal cooling lattice design procedure for injection molding applications based on expert algorithms[END_REF] optimized circular CCs in IM for dual objectives: reducing cooling time and achieving uniform cavity temperature. Four design variables, covering mold and coolant temperatures, inlet channel diameters, and spacing, were optimized. The mold's insert core was divided into hexagonal cells, each with six inlet circular CCs at cell corners and a common outlet channel. The optimization process involved solving 1D analytical equations, validated through CFD and heat transfer simulations in ANSYS Fluent and MOLDFLOW commercial software. While computationally efficient and suitable for complex geometries, this approach does not explicitly consider the 3D transient heat transfer dynamics, potentially limiting optimal solution identification. In the work of Lam et al. [START_REF] Lam | An evolutionary approach for cooling system optimization in plastic injection moulding[END_REF], a numeric tool based on MOLDFLOW coupled with the GA optimizer was developed to find the best set of variables that minimizes the temperature deviation of the cavity surface of the mold. The chromosome structure of the GA was conformed by 23 variables defining the disposition, flow rate, and coolant temperature of 3 conventional straight CCs. Because of the simplicity of the mold's and part's geometries, the boundary element method (BEM) proved to successfully accomplish the optimization task.

Nevertheless, one main limitation of the BEM is that, due to its surface-based formulation, it can lead to inaccuracies when significant heat transfer occurs within the domain's volumes or highly complex shapes are addressed. Li et al. [START_REF] Li | Topology optimization for the design of conformal cooling system in thin-wall injection molding based on bem[END_REF] also employed the BEM along with a topology optimization (TO) approach to design CCs in an IM. This scheme successfully improved efficiency and uniformity. However, the methodology does not exempt itself from relying on an initial proposal provided by the designer, and utilizes a simplified cycle-averaged approach, which may overlook temperature fluctuations ranging from 14% to 40% or more, depending on the distance between the CCs and the cavity surface [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF]. A recent study conducted by Wang et al. [START_REF] Wang | Comprehensive performance enhancement of conformal cooling process using thermal-load-based topology optimization[END_REF] introduced a comprehensive TO approach for the IM process. In their approach, the incompressible Navier-Stokes and energy conservation equations were addressed in a steady state using a pseudo-density method. The TO was performed on a prescribed 2D domain, and based on the optimized results, a 3D CCs' layout was designed. Their methodology employed a thermal-load-based approach, where a thermal heat flux proportional to the part thickness was applied to the cavity surfaces. This innovative approach allows for the consideration of parts with varying thickness. However, it is important to note that their scheme disregards the 3D heat diffusion phenomena within the part and the thermal contact resistance at the part-cavity interface, which may introduce discrepancies compared to the real heat transfer process.

In recent scientific literature [START_REF] Kanbur | Metal additive manufacturing of plastic injection molds with conformal cooling channels[END_REF][START_REF] Yun | Numerical analysis on thermo-fluid-structural performance of graded lattice channels produced by metal additive manufacturing[END_REF][START_REF] Balthazar | Conception d'un outillage de mise en oeuvre des composites à haute dynamique thermique via l'utilisation de structures lattices[END_REF], a novel methodology is proposed that involves the use of lattice structures (LSs) within the CCs of molds. The primary objective of this methodology is to enhance temperature uniformity during cooling and reduce cycle times by manipulating the porosity of the CCs using complex LSs. These structures can be customized to optimize heat flux distribution in areas of the part that are challenging to cool, such as regions with substantial thickness, while simultaneously minimizing heat flux in areas that cool rapidly, such as thin sections or corners. AM technology enables the production of these intricate and innovative geometries. In the work of Kanbur et al. [START_REF] Kanbur | Metal additive manufacturing of plastic injection molds with conformal cooling channels[END_REF], body-centered cubic (BCC) lattices were employed within tapered circular conformal CCs of a small mold made of steel using the direct metal laser sintering process. Their results showed that the introduction of six BCCs at specific locations within the tapered conformal CCs' layout yielded reductions in cooling times and temperature non-uniformity. However, it is important to note that the pressure drop increased significantly when compared to CCs without BCCs. Balthazar et al. [START_REF] Balthazar | Conception d'un outillage de mise en oeuvre des composites à haute dynamique thermique via l'utilisation de structures lattices[END_REF] arrived at similar conclusions, where BCCs and face-centered cubic (FCC) lattices were tested numerically with ANSYS Fluent on a benchmark layout reproducing IM conditions. Contrary to [START_REF] Kanbur | Design and optimization of conformal cooling channels for injection molding: a review[END_REF], where only one BCC lattice design was proposed, four different (although porosity-equivalent) geometries were addressed. Their study shows through qualitative and quantitative results that the cooling efficiency is improved by reducing the maximum, mean, and standard deviation temperatures over the working surface, potentially reducing cycle times. Further works involving LSs within CCs were addressed by Yun et al. [START_REF] Yun | Numerical analysis on thermo-fluid-structural performance of graded lattice channels produced by metal additive manufacturing[END_REF], where the porosity of FCC cells was able to change in the fluid flow direction by following three different patterns: increase-type, V-type, and W-type graded. By doing so, the heat transfer coefficient and the pressure drop can be locally modified as needed seeking for temperature homogeneity during cooling. Although the increase-type graded LS exhibited the lowest standard deviation of the working surface, it showed lower performance in terms of the average heat transfer coefficient and pressure drop compared to a constant-porosity LS. A step forward in complexity regarding AM pore structure designs for heat transfer enhancement in IM involves using triply periodic minimal surface (TPMS) structures, e.g., Schwarz-diamond and gyroid structures, as discussed in [START_REF] Oh | Adaptive conformal cooling of injection molds using additively manufactured TPMS structures[END_REF]. Their designs are generally based on mathematical cosine-sine functions allowing to generation of shells of large exchange surfaces and continuous internal channels, with the main goal of improving heat transfer efficiency while avoiding a drastic increase in pressure drop. In the aforementioned work, the authors obtained a 40% reduction in cycle time by adding the printed photo-polymer TPMS structure manufactured by the DLP technology within the CCs. While the optimization of such pore structures is in its early stages of investigation, they show great promise for improving cooling efficiency in IM processes.

Further methodologies for CCs' design can be found in the specific literature conceived for rapid tooling based on a modular-parametrical approach [START_REF] Au | A scaffolding architecture for conformal cooling design in rapid plastic injection moulding[END_REF][START_REF] Wang | Automatic design of conformal cooling circuits for rapid tooling[END_REF]. Such methodologies are based on performing an offset of the cavity surface to define a cooling design region which is then split into cells that are finally connected through lattice structures to generate the final cooling passageways. Other CCs' design methodologies exploit the offsetting concept. For instance, in the work of Zhang et al. [START_REF] Zhang | Automatic design of conformal cooling channels in injection molding tooling[END_REF], conformal surface loops are generated from the intersection of reference planes with the part surface. Then, offsets of the loops are performed which are finally connected to achieve a spiral curve that will define the CC's path. Goktas et al. [START_REF] Goktas | Cooling of plastic injection moulds using design adaptive cooling canals[END_REF] proposed a method for generating CCs in which the internal and external surfaces of the cavity are first offset by a specified distance. The edges of the offset surfaces are then used as a reference to generate equidistant paths, which are subsequently connected to produce the layout of the interior and exterior CCs.

We refer the reader to the articles [START_REF] Feng | Design and fabrication of conformal cooling channels in molds: Review and progress updates[END_REF][START_REF] Kanbur | Design and optimization of conformal cooling channels for injection molding: a review[END_REF][START_REF] Arman | A comprehensive review of injection mold cooling by using conformal cooling channels and thermally enhanced molds[END_REF] for a complete review of the aforementioned approaches. The offset-based methodologies successfully accomplished the task of assisting the mold manufacturer in the design of CCs' layouts, being able to outperform conventional straight configurations in terms of reduction of cycle time and cavity temperature homogeneity improvement. However, contrary to optimization-based design schemes [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF][START_REF] Reyes | Inverse heat transfer optimization of stamping with overmolding process involving high performance thermoplastic composites: experimental validation[END_REF][START_REF] Agazzi | Contribution à l'optimisation de la thermique des outillages dans le procédé d'injection des polymères thermoplastiques[END_REF], the thermal performance obtained by such approaches is limited, since being a purely geometrical approach, the complex physical heat transfer phenomena involved during the injection process are not taken into consideration in the design methodology.

In this work, a novel overset-FEM based optimization scheme is proposed to find a suitable arrangement of the CCs at the fixed and mobile parts of the PI. The algorithm is based on decomposing the transient heat transfer equations into three simple-to-mesh subdomains accounting for the injected part, the insert, and the CCs, and then using a high-order interpolation algorithm to obtain a global solution of the system.

The accuracy of the coupling algorithm has been already tested by the authors of this work in heat transfer problems involving both stationary [START_REF] Storti | Optimization of an internal blade cooling passage configuration using a chimera approach and parallel computing[END_REF][START_REF] Storti | A matrix-free Chimera approach based on Dirichlet-Dirichlet coupling for domain composition purposes[END_REF] and transient [START_REF] Storti | On the performance of a Chimera-FEM implementation to treat moving heat sources and moving boundaries in time-dependent problems[END_REF][START_REF] Álvarez Hostos | Solving heat conduction problems with a moving heat source in arc welding processes via an overlapping nodes scheme based on the improved element-free Galerkin method[END_REF] scenarios. The scheme allows the generation of high-fidelity structured meshes which can be anisotropically refined towards the regions where thermal gradients are higher, such as the part-mold and the CC-mold interfaces. The overset scheme is coupled to the stochastic algorithm Augmented Lagrangian Particle Swarm Optimizer (ALPSO) of the open source code pyOpt [START_REF] Perez | pyOpt: A Python-based object-oriented framework for nonlinear constrained optimization[END_REF], where 23 variables are proposed to define the layout of the CCs of the mobile and fixed parts of the mold, and an objective function term accounting for temperature homogeneity of the part and the thermal gap between maximum and minimum temperatures. The appeal of the proposed optimization technique relies on the fact that only the CC must be remeshed on each objective function evaluation. The open-source software PETSc-FEM [START_REF] Sonzogni | A parallel finite element program on a Beowulf cluster[END_REF][START_REF] Storti | PETSc-FEM, A General Purpose, Parallel, Multi-Physics FEM Program, International Center of Computational Method in Engineering (CIMEC)[END_REF], which is widely employed in the research community, is employed for both the parametric and the optimization studies that are performed throughout this communication.

This article is structured as follows. In Section 2, the study case of an actual industrial part produced by a conventional multi-component steel mold with conformal CCs is introduced. In Section 3, a comprehensive 3D numerical model is described in detail and validated, and the heat transfer of the injection process is studied. Furthermore, the thermal performance of the steel insert is compared numerically with that obtained by an identical polymer-made insert and several parametric studies are performed. Subsequently, in Section 4 a 3D optimization methodology based on overset meshes and a reduced computational model is proposed to optimize the CCs' layout in order to improve the temperature homogeneity of the part during the cooling stage. Then, in Section 5 the results of the optimization process are presented and discussed.

Finally, in Section 6, the main conclusions and future works are presented. The main limitations of employing polymer-based molds with respect to using a steel molds are also given in the aforementioned section.

Case of study

In this section, the main features regarding the geometry of the injected part under study and the multiplate injection mold are introduced. The piece considered in this work is an actual plastic part manufactured industrially by partners of this project and it is represented in Fig. 1. The bounding box which encloses the piece entirely is of dimensions 118 × 130 × 6mm 3 . Even though the geometry of the part does not present great complexity, its thickness varies markedly from 6mm at the borders, to 3.5mm towards its center (see section view of Fig. 1). Such a variation represents a thickness ratio of roughly 1.7:1 which usually leads to a non-uniform cooling of the part due to the lower rate of heat dissipation at the higher-thickness regions.

Certainly, such non-uniform cooling of the injected part eventually leads to the so undesired effects of uneven shrinkage and warpage of the manufactured piece. In addition, the studied part is only symmetric on one axis, which makes it more suitable for the application of the optimization methodology proposed in this work.

In Fig. 2 are shown the main components of the single-cavity mold used in the industry to produce the plastic part addressed in this study. As shown in the figure, the multiple-plate mold of this study is conformed mainly by the core back plate (ejector block), core plate (mold's core), and the halve insert of the mobile part, while the fixed part comprises the cavity back plate (cover block), the cavity plate (mold's cavity), and the insert cavity part. Our main objective in this work is to numerically study the feasibility, in terms of thermal performance, of replacing the conventional steel insert (SI) shown in Fig. 2d with a thermoplastic one made of polycarbonate reinforced with 10% of carbon fibers, referred to henceforth as PC10%CF. As thermoplastic materials fall into the category of plastics, which in turn fall under polymeric materials, in this paper we will refer to the PC10%CF insert with either of these terminologies, or simply PI. The temperature regulation of the plates is achieved through a series of conventional straight CCs, each with a diameter of d = 9mm, as depicted in Figs. 2b and2c. In contrast, the CCs implemented in the insert adhere to a conformal cooling concept, maintaining a consistent diameter, as depicted in Fig. 2d.

The disposition of the CCs between the fixed and mobile parts is shifted 90 • and the number of branches is different due to the geometrical constraints of the IM process, such as the positions of the ejectors and fixation holes. The channels are placed on the same plane, being its distance with the cavity surface L P -CC (see the detail view of Fig. 3) of roughly 21mm. Furthermore, in Fig. 3 it is shown a section of the entire mold assembly where can be clearly observed the location and separation of the CCs inside each component, and also the non-symmetrical disposition of the components. Because the part under study is not large enough to generate significant variations in the coolant temperature (further discussed in Section 3.4.3), a series channel arrangement is employed for the fixed and mobile parts of the steel insert, and this configuration is maintained for the plastic insert design as well. Then, the pitch distances (space between serial CCs' axes, see the detail of Fig. 3) are 38mm and 43mm for the fixed and mobile parts of the insert, respectively.

Computational model of the injection molding process

A comprehensive and accurate numerical model is introduced in this section considering the most relevant aspects of the injection molding manufacturing process. The FEM-based model comprises multiple 3D domains, multiple thermal contact resistances, and multiple transient injection cycles. The relevant process parameters are obtained from the partners of the project and described in detail. The thermal properties of the polymeric material to be used in the simulations are obtained from measurements performed on printed samples of PC10%CF.

Heat transfer governing equations

At the beginning of the injection process, the molten polymer is injected at a temperature T inj , which is superior to the polymer fusion temperature, inside the mold's cavity. Then, its heat is subtracted by conductive-convective mechanisms through the mold and the coolant fluid until reaching a desired ejection temperature T eje for the part. After the part is ejected, the mold closes, the polymer is injected, and a new cycle begins. Such a transient-cyclic thermal behavior is often modeled by averaged stationary approaches [START_REF] Xu | The design of conformal cooling channels in injection molding tooling[END_REF][START_REF] Pignon | Improvement of heat transfer analytical models for thermoplastic injection molding and comparison with experiments[END_REF]. Since temperature fluctuations around average values can be as high as 40%, complete transient heat transfer models are recommended when performing optimization [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF], and is the approach followed in this work.

For the numerical heat transfer problem, the equations to be solved are those of the 3D transient heat conduction problem with convective boundary conditions. The heat equations are solved on four domains (see Fig. 4) accounting for the injected polymer Ω 1 , the insert Ω 2 , the mold Ω 3 , and the the cover block Ω 4 .

Then, the equations to be solved are the following:

ρ 1 Cp 1 ∂T 1 ∂t = ∇ • (κ 1 ∇T 1 ) ∀T 1 ∈ Ω 1 × (0, t f ], (1)
ρ i Cp i ∂T i ∂t = ∇ • (κ i ∇T i ) ∀T i ∈ Ω i × (0, t c ], (2)
where ρ i , Cp i and κ i stands for the density, heat capacity, and the thermal conductivity of the Ω i domain, T i is the unknown temperature of each domain, and t f and t c stand for the cooling and cycle times, respectively.

It is worth noticing that Eqn. ( 1) is solved only during the cooling stage, while Eqn. ( 2) is solved during the complete cycle for the molds i = {2, 3, 4}.

The injection phase is modeled as an instantaneous non-perfect contact between the molten polymer at a temperature of T inj , and the insert at a temperature obtained from the previous cycle. Such an approach is allowed since the duration of the filling stage is negligible with respect to the other stages [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF][START_REF] Sobotka | Parametric model for the analytical determination of the solidification and cooling times of semi-crystalline polymers[END_REF]. In this context, the following initial conditions are employed: where

T 1 t = 0 + = T inj in Ω 1 , (3)
T i t = 0 + = T i t = t - c in Ω i , for i = 2, 3, 4, (4)
T i (t = 0 + ) = T i (t = t - c
) means that the temperature field T i of Ω i at the beginning of a given cycle t = 0 + is considered to be the same as the temperature field at the final time of the previous cycle t = t - c .

For the initial injection cycle, a steady-state solution is used as the initial condition in the molds, where a convective boundary condition (discussed further in the following section) is applied to the cavity's surface instead of considering the domain of the injected part.

Definition of boundary conditions

Next are described the boundary conditions applied to each domain. Due to the non-perfect contact as:

-κ1 ∂T1 ∂n = T1 -T2 T CR1-2 at Γ1-2 ∀t ∈ [0, t f ], (5)
where T CR 1-2 corresponds to the thermal contact resistance at the part-cavity (Γ 1-2 ) interface (see Fig. 4), and n stands for the outward normal of the Γ 1-2 boundary. Analogously, Eqn. (6a) defines the heat flux being absorbed by the insert through the cavity interface during cooling. Although Γ 1-2 is identical in geometry for the Ω 1 and Ω 2 domains and is therefore treated as a single entity in Eqns. ( 5) and (6a), a differentiation will be discussed in Section 3.4.1 based on the non-perfect contact assumption.

-κ2 ∂T2 ∂n =                    T2 -T1 T CR1-2 at Γ1-2 ∀t ∈ [0, t f ], (6a)
heq(T2 -Tenv), at Γ1-2 ∀t ∈ [t f , tc], (6b) T2 -T3 T CR2-3 , at Γ2-3 ∀t ∈ [0, tc], (6c)
T2 -T4 T CR2-4 , at Γ2-4 ∀t ∈ [0, tc], (6d)
After the cooling phase is finished, the part is ejected, and the Eqn. (6b) models the heat exchanged from the cavity surface to the environment during the ejection phase (t f ≤ t ≤ t c ), where T env stands for the ambient temperature, and h eq is the heat exchange coefficient which accounts for the convection and radiation. The Eqns. (6c) and (6d) accounts for the heat exchanged at the Γ 2-3 and Γ 2-4 interfaces, respectively, where T CR 2-3 and T CR 2-4 denote the thermal contact resistances at the corresponding interfaces. Furthermore, and henceforth, n should be read as a vector that is perpendicular to the boundary surface mentioned on the right-hand side of the equation, and it points away from the interior region of the domain considered on the left-hand side. Similarly to Eqns. ( 5) and ( 6), the heat fluxes at the contact boundaries of Ω 3 and Ω 4 domains are modeled as:

-κ3 ∂T3 ∂n =        T3 -T2 T CR2-3 at Γ2-3 ∀t ∈ [0, t f ], (7a)
T3 -T4 T CR3-4 at Γ3-4 ∀t ∈ [0, t f ], (7b)
and

-κ4 ∂T4 ∂n =        T4 -T2 T CR2-4 at Γ2-4 ∀t ∈ [0, t f ], (8a)
T4 -T3 T CR3-4 at Γ3-4 ∀t ∈ [0, t f ], (8b)
respectively. The heat flux exchanged between the CCs and each i-th domain (for i = 2, 3, 4) is modeled as:

-κ i ∂T i ∂n = h c (T i -T c ), at Γ i(cc) ∀t ∈ [0, t c ], (9)
being T c the temperature of the cooling fluid, h c the convective heat transfer coefficient, and Γ i(cc) denotes the boundary between each domain and the CCs, as shown in Fig. 4. A temperature uniformity assessment of the coolant fluid during the cooling stage (further discussed in Sections 3.4.3 and 6.4) shows that such variation is negligible, and allows the usage of a constant temperature of the coolant. Finally, the heat transferred from the i-th mold to the environment by convection and radiation is modeled with the following third-type condition:

-κi ∂Ti ∂n = heq(Ti -Tenv), at Γ i(out) ∀t ∈ [0, tc], (10)
where Γ i(out) is the boundary of the domain in contact with the surrounding air (see Fig. 4). It is worth clarifying that Eqn. ( 10) is valid only for Ω 3 and Ω 4 , since are the sole domains in contact with the ambient air (see Fig. 4).

Conformal grid used for FEM analysis

In Fig. 5 it is shown the conformal mesh employed for the heat transfer numerical studies. The 3D mesh is non-structured and it is conformed by a total of 6.442.206 tetrahedral linear elements, comprising a total of 1.169.447 degrees of freedom (DOF). A mesh convergence analysis was performed in order to choose the most suitable discretization parameters. The element size of the mesh is roughly 0.5mm near the part and the CCs' surfaces, where the highest thermal gradients are developed, while the mesh is coarsened towards the exterior boundaries of Ω 3 and Ω 4 domains, reaching an element size of 10mm. A detailed convergence study and a validation of the numerical model compared to analytical solutions are reported in Appendix A. It is worth mentioning that wherever a TCR is set between two domain surfaces, an algorithm was developed in the framework of this work to duplicate the nodes on each domain in order to properly model the temperature gap due to the non-perfect contact. Contrary to the work of Chen et al. [START_REF] Chen | Three-dimensional transient finite element cooling simulation for injection molding tools[END_REF], a conformal mesh approach is followed here, and it is further discussed in Section 3.4.1.

Time step adoption criterion

The sudden contact between the molten polymer and the colder cavity surface of the mold during the injection stage can give rise to numerical instabilities when modeling this phenomenon [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF][START_REF] Fachinotti | Linear tetrahedral finite elements for thermal shock problems[END_REF]. Since a refined mesh is employed in this work at the cavity-part interface along with the unconditionally stable Backward-Euler temporal integration scheme, the time step must be chosen accordingly to avoid thermal shock due to the steep thermal gradients developed in the affected region. By using tetrahedral elements and an implicit time-stepping scheme, a sufficiently high time step which accounts for the penetration depth condition must be employed, which yields [START_REF] Fachinotti | Linear tetrahedral finite elements for thermal shock problems[END_REF]:

∆t ≥ β ρ C p κ ∆x 2 , ( 11
)
where β is a constant of order 1, and ∆x is the element size at the contact region. By considering the SI properties, an element size of h = 0.5mm, and β equal to the unity for linear tetrahedral elements [START_REF] Fachinotti | Linear tetrahedral finite elements for thermal shock problems[END_REF][START_REF] Jaouen | Three-dimensional modeling by finite elements for the thermomechanical analysis of the cast parts cooling[END_REF], the minimum time step required to avoid thermal shock is ∆t ≥ 0.023s. Due to the markedly lower thermal diffusivity of the PI, the minimum required time step to avoid numerical instabilities increases to ∆t ≥ 1.162s. Regardless of the unconditionally stable temporal scheme, a compromise decision between accuracy and computational cost must be considered when adopting ∆t. Based on a numerical study performed to assess the impact of the time step adoption on the part and cavity temperatures (please see Appendix B), a ∆t = 1s was adopted for the cooling (∆t f ) and ejection (∆t e ) stages of the SI. Since different cooling and ejection times are considered for the PI throughout the parametric and optimization studies in this work, and considering that ∆t increment should be adapted on each case to the cycle time, since only integer numbers of time steps can be handled by the FEM code, such values will be reported on each particular studied case.

Thermophysical properties and injection molding parameters

With the exception of the PC10%CF material, the remaining thermophysical properties were taken from datasheets given by the mold manufacturer, which is a partner of the project, and are depicted in Table 1.

The injected material is considered to be the ABS P2MC provided by ELIX, which is an amorphous material commonly used in the manufacturing industry because of its good mechanical properties. The material properties of the steel types used in the multiplate-mold are presented in Table 1. Average properties of the injected and insert polymer materials were utilized, taking into account their working temperature range. While this assumption reduces the solution accuracy, incorporating non-linearities in the FEM studies would significantly increase the computational cost due to the large meshes, transient analysis, and multiple injection cycles involved. The main process parameters employed for the numerical simulations are depicted in Table 2. The cooling and cycle times, as well as the coolant temperature and the melted polymer injection temperature, were given by the partners of the project. The heat transfer coefficient h c between the molds and the CCs is estimated to be that according to a turbulent fluid (further discussed in Section 3.4.4), while the heat

κ[W/(mK)] ρ [kg m -3 ] Cp [J/(kgK)] α [m 2 s -1 ] Ω1 -ABS (ELIX P2MC) 0.20
T inj [ • C] h c W m 2 K t c [s] t f [s] h eq W m 2 K T env [ • C] T c [ • C] 240 10100 46 37 10 20 65
Table 2: Main parameters of the injection molding process used for the conventional steel insert provided by partners of the project.

transferred by natural convection with the ambient air and by thermal radiation is modeled through an equivalent coefficient h eq [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF], contemplating both transfer phenomena.

Handling of thermal contact resistances

When two bodies at different temperatures come into contact the temperatures at the interface of each body differ due to the existence of a heat transfer coefficient (the reciprocal of a TCR) that limits the heat flow transferred from one body to the other [START_REF] Liu | Evaluation of heat transfer coefficient between polymer and cavity wall for improving cooling and crystallinity results in injection molding simulation[END_REF]. This coefficient depends on factors such as the surface rugosity of the tool, the injection pressure, and the wettability of the polymer, among others [START_REF] Somé | Modeling of the thermal contact resistance time evolution at polymer-mold interface during injection molding: Effect of polymers' solidification[END_REF].

Since non-perfect contacts are considered at the mold-mold and part-mold interfaces in this work, two main approaches can be followed in the FEM context to cope with the different temperatures of each contact surface: the conformal and the non-conformal mesh approaches [START_REF] Chen | Three-dimensional transient finite element cooling simulation for injection molding tools[END_REF] (see Fig. 6). For complex geometries, the non-conformal mesh scheme followed by Chen et al. [START_REF] Chen | Three-dimensional transient finite element cooling simulation for injection molding tools[END_REF] (see Fig. 6a) for the simulation of a complete IM process is more practical since each domain can be meshed independently. However, a projection is required due to the non-matching nodes, which induces an asymmetry of the system matrix that can affect the convergence of the iterative solver. To facilitate the conformal scheme shown in Fig. 6b, an algorithm was devised within the scope of this study to automatically duplicate nodes at the contact surfaces. Since the domain boundaries are duplicated because of the addition of new nodes by the conformal approach proposed here (see Fig. 6b), we will use henceforth the notation Γ i-j to refer to the boundary of Ω i in contact with Ω j (red nodes of Fig. 6b). Analogously, Γ j-i will be used to refer to the boundary of Ω j in contact with Ω i (green nodes of Fig. 6b).

without TCR with TCR

Ω j Ω i Ω j Ω i (a)
Non-conformal approach (Chen et al. [START_REF] Chen | Three-dimensional transient finite element cooling simulation for injection molding tools[END_REF]).

without TCR

with TCR To describe how the different TCRs are handled, let us consider the following algebraic form of Eqn. [START_REF] Davoudinejad | Experimental investigation and thermomechanical modelling for tool life evaluation of photopolymer additively manufactured mould inserts in different injection moulding conditions[END_REF] after a standard FEM Galerkin formulation is applied:

Ω i Ω j Γ i-j Ω i Ω j Γ i-j Γ j-i (
C Ṫ + KT = Q, ( 12
)
where C, K and Q are the global capacitance matrix, conduction matrix, and heat load vector, respectively, which are defined as:

C = Ω ρCpN T N dΩ, K = Ω κB T B dΩ, Q = Γ q N T dΓ, (13)
being B ij = ∂N i /∂x j , such that ∇T ≈ BT, and the temperature field is approximated as:

T (x, t) ≈ Nj(x)Tj(t) = N(x)T(t) for all x ∈ Ω, ( 14
)
where N j is the shape function associated to the node j of the discretized finite element space of Ω, and T j the unknown temperature of the node. Now, by using the unconditionally stable implicit backward Euler scheme to integrate the transient heat transfer problem of Eq. ( 12), the algebraic problem yields:

C T t+∆t -T t ∆t + KT t+∆t = Q t+∆t , ( 15
)
where ∆t is the adopted time step. In a first instance assume that Ω is isolated, i.e. q = 0 at Γ. Then, by rearranging the terms of Eqn. ( 15) we obtain:

1 ∆t C + K A Ω T t+∆t = 1 ∆t C T t rhs Ω , (16)
where A Ω refers to the thermal system matrix, and rhs Ω the right-hand side of the algebraic equation.

Then, extending Eqn. [START_REF] Rahmati | Rapid tooling analysis of Stereolithography injection mould tooling[END_REF] for the multiple domains involved in the IM process, the following system is obtained:

      AΩ 1 0 0 0 0 AΩ 2 0 0 0 0 AΩ 3 0 0 0 0 AΩ 4             TΩ 1 TΩ 2 TΩ 3 TΩ 4       t+∆t =       rhsΩ 1 rhsΩ 2 rhsΩ 3 rhsΩ 4       , ( 17
)
and consequently, by considering the heat transfer through the surfaces, the algebraic system to be solved yields:

      A ′ Ω 1 -AΩ 1 Ω 2 0 0 -AΩ 2 Ω 1 A ′ Ω 2 -AΩ 2 Ω 3 -AΩ 2 Ω 4 0 -AΩ 3 Ω 2 A ′ Ω 3 -AΩ 3 Ω 4 0 -AΩ 4 Ω 2 -AΩ 4 Ω 3 A ′ Ω 4             TΩ 1 TΩ 2 TΩ 3 TΩ 4       t+∆t =       rhsΩ 1 rhs ′ Ω 2 rhs ′ Ω 3 rhs ′ Ω 4       , ( 18
)
where A ΩiΩj stands for coefficients matrix added due to the heat flux at the contact of Γ i -Γ j , and can be read as:

AΩ i Ω j = Γ i-j 1 T CRi-j N T NdΓ. ( 19
)
Furthermore, due to the coupling between domains by the TCR condition, and the heat flux transfer to the CCs and/or the ambient (when apply), a mass matrix M Ωi is added to the A Ωi matrix of Eqn. [START_REF] Bagalkot | Categorization of failures in polymer rapid tools used for injection molding[END_REF], becoming A ′ Ωi in Eqn. [START_REF] Kovács | Thermal simulations and measurements for rapid tool inserts in injection molding applications[END_REF]. The mass matrix obtained for each domain is defined as:

M Ω1 = Γ1-2 1 T CR 1-2 N T NdΓ, (20)
M Ω2 = Γ 2(CC) h c N T NdΓ + F Ω4 = - Γ 4(CC) h c T c N T dΓ - Γ 4(ext) h eq T env N T dΓ, (26)
which, after adding them to the right-hand side of the equation system, the modified rhs ′ Ωi vector is obtained for each domain. Then, the systems are solved in a monolithic manner in PETSc-FEM [START_REF] Storti | PETSc-FEM, A General Purpose, Parallel, Multi-Physics FEM Program, International Center of Computational Method in Engineering (CIMEC)[END_REF].

Regarding the thermal contact resistance between the ABS polymer and the steel insert T CR 1-2 , a constant value of 10 -3 m 2 K W -1 is considered in this communication [START_REF] Pignon | Improvement of heat transfer analytical models for thermoplastic injection molding and comparison with experiments[END_REF]. For the contact between the steel molds, a constant value of 3.33 × 10 -5 m 2 K W -1 [START_REF] Chen | Three-dimensional transient finite element cooling simulation for injection molding tools[END_REF] is employed. Given the lack of specific information in the literature on TCR for newly developed polymer matrix composites [START_REF] Zhou | Fundamental and estimation of thermal contact resistance between polymer matrix composites: A review[END_REF][START_REF] Le Mouellic | Development of an original overmoulding device to analyse heat transfer at polymer/polymer interface during overmoulding[END_REF], TCR values of the plastic insert are approximated as follows: the TCR for the contact between the PC10%CF insert and the steel mold is taken to be 6 × 10 -3 m 2 K W -1 [START_REF] Fuller | Thermal contact conductance of metal/polymer joints: an analytical and experimental investigation[END_REF], and the polymer-polymer interface is modeled considering a value of 10 -2 m 2 K W -1 [START_REF] Le Mouellic | Development of an original overmoulding device to analyse heat transfer at polymer/polymer interface during overmoulding[END_REF].

Measurements of thermal properties of the FDM insert material

The thermal conductivity of the 3D printed composite material PC10%CF was measured by the hot guarded plate methodology [START_REF] Goff | Thermal conductivity of an injected polymer and short glass fibers composite part: measurement and model[END_REF] whose working principle is mainly based on a one-dimensional analysis of the Fourier law in stationary regime. To this end, a homemade device located at the laboratories of Laboratoire de Thermique et Énergie de Polytech, Nantes Université, was employed. The dimensions of the printed samples were 15 × 15 × 3mm 3 and the exterior dimensions of the guard were 40 × 40 × 3mm 3 , and its inner dimension 20 × 20 × 3mm 3 . Both the samples and the guards were produced by a Lynxter S600D printer device, and their infill was 100%. The thermal conductivity measurements were performed at temperatures from 40 • C to 125 • C. However, no temperature dependency of the material thermal conductivity was observed during measurements. The temperatures of the hot and cold sides of the measurement device were controlled by cooling liquids whose temperatures were set by two thermoregulators LAUDA RP 855. The specific heat of the PC 10%CF model material was determined at the laboratory as well, with a Q200 Differential Scanning Calorimetry (DSC) device. The specific heat of the PC 10%CF material was analyzed in the range 20 • C -285 • C with single heating. The heating rate was 10 • C min -1 and the mass of the sample was 35mg.

From the DSC measurements, the T g of the sample was identified to be at 142.2 • C. The value of the heat capacity presented in Table 1 corresponds to a temperature of 100 • C, and is considered to be constant in this work.

Temperature uniformity of the coolant

Initially, the coolant is assumed to be water at 65 • C, and its flow rate is given by the thermoregulator operational condition, which is considered to be in this work of VF = 9lmin -1 on the fixed and mobile parts. By performing a heat balance assuming that all the heat subtracted from the part during cooling is transferred to the CC, i.e., by neglecting the heat loss through the external mold boundary [START_REF] Xu | The design of conformal cooling channels in injection molding tooling[END_REF][START_REF] Park | Optimal cooling system design for the injection molding process[END_REF], the heat balance within the mold yields:

QP + QCC = 0, ( 27
)
where QP is the rate of heat flow from the melted part estimated for an amorphous polymer as [START_REF] Rao | Optimization of cooling systems in injection molds by an easily applicable analytical model[END_REF]:

QP = - 1 t f VP ρP Cp P Tinj -T eje , (28)
being v p and ρ p the volume and density of the part, and QCC being the rate of heat flow exchanged with the coolant computed as:

QCC = 10 -3 60 VF ρF CpF (Tc o -Tc i ) ∆Tc i-o , ( 29
)
where ∆T ci-o is the temperature difference between the inlet temperature T ci and the outlet temperature T co , while ρ F and Cp F are the density and the heat capacity of the fluid. From such estimation it is obtained that the temperature variation of the coolant is negligible, (∆T ci-o ≈ 0.3 • C), allowing the simplification of the temperature variation inside the CCs to a constant value, which is a common practice in the IM process modelling [START_REF] Feng | Design and fabrication of conformal cooling channels in molds: Review and progress updates[END_REF][START_REF] Chen | Three-dimensional transient finite element cooling simulation for injection molding tools[END_REF][START_REF] Torres-Alba | A hybrid cooling model based on the use of newly designed fluted conformal cooling channels and fastcool inserts for green molds[END_REF]. Even though the aforementioned simplification is employed throughout this work, it should not be considered as a limitation of the optimization methodology (later depicted in Section 5), since the scheme could be directly extended to solve diffusion-advection inside the CCs to account for large variations of ∆T ci-o , or, furthermore, to be coupled to the Navier-Stokes module of PETSc-FEM for a complete conjugated heat transfer approach. This topic will be further discussed in the section 6.4.

Heat transfer coefficient in the CCs

The convective heat transfer coefficient between the bulk of the coolant fluid inside the channels and the mold surface is estimated as follows:

hc = κF Nu d , (30)
where κ F is the thermal conductivity of the coolant, d is the hydraulic diameter of the passage, and Nu stands for the Nusselt number which is obtained according to the Dittus-Boelter correlation for turbulent flows, which states:

Nu = 0.023 Re 0.8 Pr 0.4 , (31)
where Pr and Re stand for Prandtl and Reynolds numbers, respectively, which are defined as:

Re = ρF d vF µF , Pr = CpF µF κF , (32)
where ρ F , µ F , and v F are the density, dynamic viscosity, and velocity of the fluid, respectively, and Cp F is the heat capacity of the coolant at constant pressure. Consequently, the convective heat transfer coefficient is estimated to be h c = 10100W/(m 2 K), which is in agreement with reported values usually employed in injection molding processes [START_REF] Sobotka | Optimisation thermique des outillages d'injection thermoplastiques[END_REF] (10000W/(m 2 K) -15000W/(m 2 K)). It is worth clarifying that the Dittus-Boelter correlation is only valid within specific ranges of the Re, Pr, and length-to-diameter ratio values [START_REF] Storti | Optimization of an internal blade cooling passage configuration using a chimera approach and parallel computing[END_REF], which are all well-suited in the present communication.

Thermal performances of the steel and polymeric inserts

In this section, 3D FEM analyses are performed to assess and compare the thermal performance of the insert previously depicted by considering the conventional steel material, and a polymeric material. The suitability of the numerical model is assessed by an experimental investigation using a half-plastic half-steel insert. Furthermore, one-dimensional FEM simulations are performed to study the sensibility of the process cycle time to parameters such as the thermal conductivity of the polymer, distance from the CC to the injected part, and the thermal contact resistance at the cavity-mold interface.

Thermal performances of the steel and polymeric inserts

The results of the temperature obtained by the 3D heat transfer solutions at the part and the insert are presented in Fig. 7. The governing equations were solved with the advection-diffusion module of PETSc-FEM by means of the Conjugated Gradient (CG) iterative solver. The time increment for the cooling and the ejection stages of the PI were ∆t f = 3.7s and ∆t e = 3s, respectively. The mean temperatures of a given domain Ω i , and the mean temperature at its surface Γ i at a time t, are computed from the FEM solution as follows:

T Ωi (t) = 1 Ω i Ωi T i (t)dΩ i ≈ 1 Ω i n el j=1 V j i 4 4 k=1 T k i (t) , (33)
T Γi (t) = 1 Γ i Γi T i (t)dΓ i ≈ 1 Γ i n el j=1 A j i 3 3 k=1 T k i (t) , (34)
where n el are the number of elements of a given i-th domain, and V j i and A j i correspond to the volume and surface of a tetrahedral/triangular j-th element, respectively.

It can be observed from the results of Fig. 7 that the thermal performance of the PI greatly differs from that of the SI for the same injection process parameters. Certainly, the difference between their thermal diffusivity plays a critical role in cooling the injected part. It is observed from Figs. 7a and 7b that the periodic steady-state (PSS) is quickly reached after the fourth cycle when using the steel insert, achieving a surface temperature of the part of ≈ 74.9 • C at the end of the cooling stage. This value aligns with the demolding temperature requirement specified in the datasheet provided by the mold manufacturer, who is a collaborating partner in the SAMFAST project, indicating a demolding temperature of 75 • C for ABS.

The result of the FEM model is also consistent with results obtained by an analytical approach depicted in Appendix A. Furthermore, the mean temperature of the part volume is ≈ 95.7 • C after cooling, which is below the glass transition temperature of the ABS material (T g = 105 • C), recommended condition to eject the part.

Nevertheless, as can be seen from Figs. 7c and7d, the scenario for the thermoplastic insert is completely different, since more than 20 cycles would be required to reach a PSS of the process. It is important to note that, in this study, the attainment of the PSS is determined by the condition that the mean temperature difference on the cavity surface between two successive injection cycles is sufficiently small, with a criterion set at a temperature difference of 0. In Fig. 8 are shown the temperature distribution at the insert surface (Γ 2-3 ), and at the mold surface (Γ 3-2 ), which are in contact with each other, at the final of the ejection time once the PSS is reached.

Due to the non-perfect consideration between both surfaces, such temperatures are not equal. On the one hand, it is observed from Figs. 8a and 8b that for the SI, the temperature distributions of both surfaces are similar. On the other hand, for the PI, a marked difference in the temperature distribution is observed between the surfaces shown in Figs. 8c and 8d due to the higher value of the TCR and the great contrast between the thermal properties of the materials in contact. Given the relatively stable nature over time of the temperature field obtained at the mold and block surfaces which are in contact with the PI, i.e. at Γ 3-2 ∪ Γ 4-2 , such a temperature field will serve as a suitable basis for defining the BCs in the optimization task discussed in Section 5.

Experimental validation of the numerical model

The numerical model demonstrated its suitability when applied to the conventional SI case, and the results were compared to the corresponding datasheet. However, before proceeding to the subsequent parametric and optimization investigations that constitute the primary focus of this work, it is essential to conduct further validation of the numerical model for the PC10%CF case. In this regard, the mobile part of the PI insert with the original CC's layout was manufactured using FDM. An experimental validation involving a half-plastic, half-steel insert configuration is presented next.

Setup of the experience

Experimental essays were performed on an IM machine model ZERIES ZE1900, shown in for consistent results. To enable an accurate comparison, two simulation scenarios were considered: one involving the presence of four steel ejectors within the plastic half (depicted in Fig. 9b), denoted as FEM (WE), and another scenario where the ejectors were disregarded, denoted as FEM (WOE). For the numerical simulations, ∆t f = 2s was adopted.

Following the cooling phase, thermal images were captured using a FLUKE TiS75 thermal camera to assess the temperature distribution of the mold and the injected part. To ensure precise temperature readings, the emissivity value of the infrared (IR) camera was carefully configured. In order to determine the accurate emissivity and reflectivity values of PC10%CF and ABS samples, a Bruker 80V spectrometer device was utilized at the wavelength range corresponding to the IR camera. The detailed measurements conducted for emissivity determination can be found in Appendix C. Based on the acquired results, emissivity values of 0.916 were assigned to the PC10%CF material and 0.929 to the ABS material, respectively.

Numerical and experimental comparison

Fig. 10a displays the surface of the injected part immediately after the mold is opened, while Fig. 10b illustrates the surface cavity right after the ejection of the part from the mold. In both cases, a compre- hensive comparison is presented, featuring the actual photograph, the temperature field obtained through the utilization of an IR camera, the postprocessed IR image, and the numerical solution, arranged from left to right. Notably, the IR imaging was solely performed on the PI side, and it is important to acknowledge that due to practical constraints, perpendicular imaging was not feasible, potentially introducing slight uncertainties in the obtained temperature fields.

From a qualitative perspective, it is observed in Fig. 10 a remarkable agreement between the temperature fields obtained at the surface of the injected part, which is in contact with the SI during the cooling process, and the cavity surface of the PI, when compared to the FEM (WE) simulation conducted using the numerical model presented in this study. The images captured by the IR camera were subjected to postprocessing in Paraview to ensure consistency with the colorbars employed for the simulation postprocess, and to facilitate the line plots performed in this section. It can be observed that the locations of maximum temperature at both the part and cavity surfaces observed in the IR camera images are in good agreement with those obtained by the FEM solution. A quantitative analysis of the results, as illustrated in Table 3, confirms a high degree of agreement between the numerical and experimental data. The T max Γ1-2 and T Γ1-2 values exhibit a difference of approximately 2 -3 • C between the FEM (WE) solution and the temperatures obtained from the IR camera measurements. Furthermore, an excellent match is observed in the temperature distribution at the surface of the injected part, with a difference of σ 1-2 below 1 Part surface (steel side) Cavity surface (plastic side) σ 1-2 values correspond to the specific region highlighted in magenta in Fig. 10. However, for the temperature distribution at the cavity surface of the PI, the differences slightly increase to approximately 2

T max Γ 1-2 [ • C] T Γ 1-2 [ • C] σΓ 1-2 [ • C] T max Γ 2-1 [ • C] T Γ 2-1 [ • C] σΓ 2-1 [ • C] IR camera
• C -3 • C,
indicating a slight overestimation of the numerical simulation. This discrepancy can be attributed to several factors, including the assumption of constant thermophysical properties of the polymers with respect to temperature, the accuracy of the TCR value, and the quality of the printed PI material.

The FEM (WOE) case, as indicated in Table 3, demonstrates a very good performance with respect to experimental measurements as well. The T Γ2-1 value is nearly 3 • C higher compared to the FEM (WE) case and exceeds the experimental results by ≈ 5 • C. This disparity is expected since the additional heat dissipated by the steel ejectors is disregarded in the FEM (WOE) scenario. However, the remaining reported values exhibit a high level of agreement with the experimental data, particularly the standard temperature deviations at the surfaces of the part and insert. These results hold particular importance as they provide evidence that the numerical model exhibits minimal error when neglecting the ejectors. Further comparisons are shown in Fig. 11 where the temperatures are plotted over vertical (A-A) and horizontal (B-B) lines placed at the geometrical center of the part and cavity surfaces. Good agreement of the temperature distributions is observed as well. Based on the above, the numerical model is considered accurate enough for conducting the parametric and optimization studies described in subsequent sections of this article.

Sensitivity and order-of-magnitude analysis

The optimization of heat transfer in an IM process involves modifying specific physical and geometrical parameters to facilitate efficient heat dissipation from the injected polymer to the coolant. Despite the extensive literature on parameter selection for steel molds, there is limited research focusing on a thorough comprehension and appropriate application of these parameters for PIs [START_REF] Schuh | Towards Temperature Control Measures for Polymer Additive Injection Molds[END_REF].

By following the electrical analogy, i.e., considering the overall heat transfer on a composite wall as the current passing through several conduction and convection serial resistances [START_REF] Rao | Design formulas for plastics engineers[END_REF], and considering a 525 symmetric 1D case as shown in Fig. 12a, the heat flux density ( Q/S) transferred from the center of the part to the CCs can be estimated as: where the numerator of the right-hand side stands for the temperature difference between the center of the part T P and the coolant temperature T c , and the denominator is the global thermal resistance which accounts for the thermal resistances of the injected polymer R P of thickness e, the contact part-mold T CR 1-2 , the insert R M , and the convective resistance R C .

d Q dS = TP -Tc e/2 κ1 R P +T CR1-2 + LP -CC κ2 R M + 1 hc R C , (35)
From the order-of-magnitude estimation of the table of Fig. 12b it is observed that for a SI, the resistances leading the heat transfer are those represented by the injected part and the TCR. On the other hand, when employing a PI, the mold resistance R M becomes also crucial in the heat transfer process, which is roughly two orders of magnitude lower than that of the SI. Hence, it is of interest to study how the TCR and the R M impact on the t c for the polymer-based mold. The modification of parameters involving substantial reductions in coolant temperature and extensions of the cycle time is common in PIs for achieving desired ejection temperatures [START_REF] Davoudinejad | Experimental investigation and thermomechanical modelling for tool life evaluation of photopolymer additively manufactured mould inserts in different injection moulding conditions[END_REF][START_REF] Whatcott | Effects of Conformal Cooling Channels on Additively Manufactured Injection Molding Tooling[END_REF][START_REF] Bogaerts | Influence of thermo-mechanical loads on the lifetime of plastic inserts for injection moulds produced via additive manufacturing[END_REF][START_REF] Kovács | Thermal simulations and measurements for rapid tool inserts in injection molding applications[END_REF]. In light of the detrimental effects of high mold temperatures and hot spots on insert durability [START_REF] Bogaerts | Influence of thermo-mechanical loads on the lifetime of plastic inserts for injection moulds produced via additive manufacturing[END_REF][START_REF] Bagalkot | Categorization of failures in polymer rapid tools used for injection molding[END_REF], a coolant temperature of T c = 10 • C is adopted for the numerical simulations henceforth. .

Influence of the insert thermal resistance on the cycle time

To assess the sensitivity of t c with respect to R M , transient heat transfer problems are solved involving the Eqns. ( 1) and ( 2) on two 1D computational domains accounting for the injected polymer and the insert. The sensitivity analysis is performed by increasing the κ 2 from 0.1W/(mK) to 1.2W/(mK) in steps of 0.1W/(mK). Furthermore, two distances of the CCs are considered: L P -CC = 21mm, which is the original disposition, and L P -CC = 9mm, which is the minimum distance recommended in the literature (L P -CC ≥ d [START_REF] Kanbur | Design and optimization of conformal cooling channels for injection molding: a review[END_REF]). The remaining thermophysical properties remain unchanged (see Table 1). Then, the parametric problem consists in finding the t c for a given value of the R M , i.e., for a given set of L P -CC and κ 2 values, which minimizes the following function:

f (tc, RM ) = T ref Ω 1 -T Ω 1 (tc, RM ) 2 , ( 36
)
where T ref Ω1 stands for the mean temperature of the injected part obtained by the SI at the end of the cooling stage, being 95.7 • C (see Section 4.1), and T Ω1 (t c , R M ) is the mean temperature of Ω 1 for the study configuration computed in the FEM 1D mesh as:

T Ω 1 (tc) = 1 Ω1 Ω 1 T1(tc)dΩ1 ≈ 1 Ω1 n el j=1 L j 1 2 2 k=1 T k 1 (tc) , (37)
where L j i is the length of the j-th element of the mesh. The criterion that considers the mean temperature of the part as ejection temperature was employed in this analysis following the work [START_REF] Pignon | Improvement of heat transfer analytical models for thermoplastic injection molding and comparison with experiments[END_REF]. The 1D mesh employed for the simulations comprised a total of 964 nodes refined toward the part-insert contact point.

The time step was defined as 1s for both the cooling and ejection stages. Since a smaller ∆x was employed for the 1D mesh compared to the 3D one, the usage of ∆t = 1s was also possible for the PI. The results are presented in Fig. 13 where power-law functions were employed to fit the points. Results of Fig. 13 show that the t c of the injection process is very sensitive to both the conductivity value of the thermoplastic insert and the distance between the injected part and the CC. For instance, for an injected part of thickness properties, outlined in Table 1.

4.5mm on a PI with κ = 0.1W/(mK), t c can take values of 22.7 and 8.2 times higher than that employed by the actual SI with the reference configuration of the CC (L P -CC = 21mm) and the closer configuration (L P -CC = 9mm), respectively. Such values of t c are reduced to 6.4 and 2.4 times higher than the time required by the steel insert when κ = 0.33W/(mK) is considered. These findings carry significant implications, highlighting the need for extensive research to investigate high-performance polymers that exhibit thermal conductivities surpassing 1.0W/(mK). This is essential to ensure that PIs can maintain their competitiveness in terms of productivity when compared to SIs. In fact, thermal conductivity enhancement of polymers for heat exchanger applications is the main concern of recent research studies [START_REF] Hussain | Review of polymers for heat exchanger applications: Factors concerning thermal conductivity[END_REF][START_REF] Ahmadi | Performance Evaluation of hi-k Lung-inspired 3D-printed Polymer Heat Exchangers[END_REF].

Sensibility of the cycle time to the TCR

Further parametric studies were performed to assess the sensibility of the TCR at the part-insert interface on the cycle time. Numerical studies were performed for the PI and SI to compare such an influence in both types of insert materials. Then, the objective of the parametric study is to minimize the function of Eqn. [START_REF] Storti | Optimization of an internal blade cooling passage configuration using a chimera approach and parallel computing[END_REF], however, the mean temperature of the injected part T Ω1 (t c , R M ) now becomes T Ω1 (t c , T CR), being a function of the cycle time and the TCR value, solely. The thermophysical properties of the part and insert remain unchanged, and are those reported in Table 1. The distance between the CCs and the part surface corresponds to that of the original disposition (L P -CC = 21mm). Fig. 14 presents the dependence of cycle times required to cool down the part to the ejection temperature on TCR and part thickness. It can be observed that both the SI and PI cycle times are significantly influenced by high TCR values, emphasizing the importance of considering it for numerical simulations. However, their behavior differs, as the TCR of the SI can be disregarded only for values lower than 10 -4 m 2 K W -1 without significant error, consistent with the literature [START_REF] Pignon | Improvement of heat transfer analytical models for thermoplastic injection molding and comparison with experiments[END_REF]. In contrast, the TCR of the PI can be disregarded for values lower than 10 -3 m 2 K W -1 as the mold resistance starts to dominate heat transfer phenomena.

Moreover, the results of Fig. 14 support the choice of the TCR value for the polymer-polymer contact in Section 3.4.1, since values exceeding 10 -2 m 2 K W -1 start to lose physical significance for a solid-solid contact.

3D parametric study of the CC position

In light of a significant literature gap on the accurate assessment of thermal performance for PI with CCs, this study conducts a parametric investigation to examine how the part surface temperature is impacted by the proximity of CCs. Three discrete values, specifically 21mm (case 1), 15mm (case 2), and 9mm (case 3), are considered. 3D simulations with multiple injection cycles are executed for each scenario until reaching the PSS. To enable an equitable comparison between the PI and SI, this parametric study aims to identify the cycle time required to achieve a surface mean temperature similar to that of the SI (see Section 4.1).

To this end, a mean temperature range of [74 • C -76 • C] is utilized as a criterion. In addition, the ejection time is extended to up to 60s. A time of 5s was used for ∆t f and ∆t e in the FEM simulations.

As a result, the cooling times needed to reach the ejection temperature were 360s, 310s, and 240s for the respective cases. It can be seen from Fig. 15 how the homogeneity of the temperature field obtained at the end of the cooling stage decreases as the CCs get closer to the cavity surface. This can also be noticed in Fig. 16, where the temperature of the part surface is plotted in steps of 6 • C. From such results, it was observed that 56.95%, 48.84% and 52.25% of the part surface temperature was between the range of [65

• C -83 • C]
at the end of the cooling time for the three cases, respectively. On the other hand, such a value was 100% for the SI.

Further results are also depicted in Table 4, where the temperature homogeneity at the part surface is assessed as follows [START_REF] Tang | A Computer-Aided Optimization Approach for the Design of Injection Mold Cooling Systems[END_REF]: It is worth noting that by normalizing Eqn. [START_REF] Storti | On the performance of a Chimera-FEM implementation to treat moving heat sources and moving boundaries in time-dependent problems[END_REF] with respect to the part surface area and taking the square root, the standard deviation of temperature, denoted as σ 1-2 , can be obtained. The values of T max Γ1-2 and T min Γ1-2 in Table 4 represent the maximum and minimum temperature of the part surface at the end of cooling. It can be noticed from Table 4 that, despite the CCs' position, the temperature homogeneity of the part is severely affected by the PI with respect to the SI. It can also be observed that T max Γ1-2 increases as L P -CC is reduced, increasing also the temperature gap ∆T Γ1-2 between maximum and minimum temperatures at the part surface. These elevated temperatures are probably the cause of hot spots in the insert, identified as one of the main factors for its premature failure. Moreover, the results unveil contrasting temperature gradients within the part thickness during the cooling phase. Notably, when employing the conventional SI, a considerable temperature difference of almost 21 • C between the mean surface and volume temperatures is observed. Conversely, the utilization of the PI yields significantly diminished temperature gradients ranging from 2 • C to 5 • C. This dissimilarity arises from the PI's lower heat diffusivity, which restrains the rapid dissipation of heat from the part, thereby promoting a more uniform cooling process which can eventually lead to warpage reduction. In the next section, an optimization methodology is proposed to improve the temperature homogeneity of the part surface while reducing its temperature gradients.

T dev Γ 1-2 = Γ 1-2 T (t f ) -T (t f ) 2 dΓ. ( 38
)

Optimization methodology

Achieving temperature homogeneity in IM while maintaining a reasonable cycle time is a challenging task. The use of PI exacerbates this difficulty due to the high temperature deviation at the part surface.

To address this issue, we present a novel optimization methodology in this paper that leverages the overset-FEM technique and the ALPSO algorithm. Our methodology incorporates the flexibility of FDM technology, allowing the CCs to bend as needed. 625

T dev Γ 1-2 [ • C 2 m 2 ] 0.

Consideration of a reduced computational domain

The time-independent temperature behavior at the Γ 2-3 ∪ Γ 2-4 interfaces was established in Section 4.1. Since the optimization procedure involves changing the geometry of the CCs, it is infeasible to solve the full 3D model for each OFE. Instead, we employ a reduced computational model that considers only Ω 1 and Ω 2 , using the full model results to establish proper boundary conditions at the insert. More precisely, the transient heat transfer Eqns. ( 1) and ( 2) are to be solved in the reduced computational model for the injected polymer and the PI only, with the initial conditions of Eqns. ( 3) and (4), and BCs of Eqns. ( 5) and [START_REF] Kanbur | Metal additive manufacturing of plastic injection molds with conformal cooling channels[END_REF]. Notably, instead of solving for T 3 and T 4 , the contacts between the PI with the molds and blocks of Eqns. (6c) and (6d) are modeled here as:

-κ2 ∂T2 ∂n =          T2 -T f ull 3-2 T CR2-3 , at Γ2-3 ∀t ∈ [0, tc], (39a)
T2 -T f ull 4-2 T CR2-4 , at Γ2-4 ∀t ∈ [0, tc], (39b)
where T f ull 3-2 and T f ull 4-2 are the temperature fields fixed in time obtained from the full 3D computational model at Γ 3-2 and Γ 4-2 , respectively, at the end of cooling. A linear interpolation is performed to transfer the temperature field from the unstructured mesh employed by the full 3D model to a structured mesh used by the reduced computational model employed in the overset-based optimization scheme described next.

Overset mesh FEM methodology

To explain the overset scheme followed in the proposed optimization process, a two-dimensional scheme of the insert (Ω 2 ) and the injected polymer (Ω 1 ) is considered, as shown in Fig. 17, where the transient diffusion heat transfer equation needs to be solved. Let us assume that the reduced computational domain, i.e. Ω 1 ∪ Ω 2 , is described by the overlapping of the following domains (see Fig. 17a):

• Ω CC accounting for the CCs' description,

• Ω M comprising the entire insert domain,

• Ω P describing the surface cavity and the injected polymer.

A mesh is required to discretize each domain for the purpose of simulating the system, with thermophysical properties ρ, C p , κ assigned to each mesh element. In regions where the overset grid overlaps with Ω 1 , the mesh elements inherit properties of the melted polymer, while in zones where the grid overlaps with Ω 2 , the mesh elements have the insert properties. The domain Ω P illustrated in Fig. 17 can be decomposed into two disjoint subdomains, Ω P = Ω P 1 ∪ Ω P 2 . Then, the system of equations to be solved yields:

     A ′ Ω P 1 -AΩ P 1 Ω P 2 0 0 -AΩ P 2 Ω P 1 A ′ Ω P 2 0 0 0 0 A ′ Ω M 0 0 0 0 A ′ Ω CC           TΩ P 1 TΩ P 2 TΩ M TΩ CC      t+∆t =      rhsΩ P 1 rhsΩ P 2 rhs ′ Ω M rhs ′ Ω CC      , (40)
where it can be seen that the systems Ω P 1 and Ω P 2 are linked through the T CR condition at their interfaces, however, they are not related to the mold domain Ω M . In addition, Ω M and Ω CC domains are also unrelated to each other in Eqn. [START_REF] Perez | pyOpt: A Python-based object-oriented framework for nonlinear constrained optimization[END_REF]. In order to obtain a global solution of the heat transfer equation within the insert, the part, and the CCs, a coupling strategy must be performed to transfer information between the overlapped domains. To this end, the accurate high-order coupling scheme of Storti et al. [START_REF] Storti | Optimization of an internal blade cooling passage configuration using a chimera approach and parallel computing[END_REF][START_REF] Storti | A matrix-free Chimera approach based on Dirichlet-Dirichlet coupling for domain composition purposes[END_REF] is employed in this work, which is briefly described next for completeness of the article.

Nodes grouping

After each domain has meshed, the nodes of each domain are grouped in sets according to their locations and/or functionality. Therefore, considering a given domain ξ, three main nodes sets can be distinguished:

• set B for BC's: boundary nodes of ξ which are coincident with the boundary nodes of Ω 1,2 , i.e., B ξ = {∂Ω ξ ∩ ∂Ω 1,2 }, which are used to define the boundary conditions of the thermal problem,

• set I for interpolation: nodes located at the boundary of ξ that do not belong to Ω 1,2 boundaries, i.e.,

I ξ = {∂Ω ξ \ ∂Ω 1,2
}, which are used to interpolate values between subdomains.

• set Z of internal nodes: nodes of ξ which are neither part of set B nor of set I, i.e., Z ξ = {ξ \{B ξ ∪I ξ }},

29 Ω Ω p(Ω 2 ) Γ CC Γ m Ω CC Ω P =Ω P1 ∪Ω P2 Ω M Ω Ω p(Ω 2 ) Γ CC Γ m Ω 1 Ω 2 Ω 2

Domain decomposition

(a) Domains considered by the reduced model (left) and the domains obtained after the domain decomposition for the overset-FEM scheme. A scheme of the defined boundary nodes and the interpolation nodes of each domain can be observed in Fig. 17b, where the subdomain meshes are schematized. The remaining nodes are assumed to be those from the Z group, not plotted in Fig. 17 for clarity.

Ω M Ω CC Ω P(Ω 2 ) Ω P(Ω 1 ) Γ CC Γ 1-2 Ω 1 Ω 2 Γ M Interpolation nodes (I) of
for domains Ω M and Ω CC , with their corresponding node-splitting, then the equation yields:

                 A ′ BB Ω M A ′ BZ Ω M A ′ BI Ω M 0 0 0 A ′ ZB Ω M A ′ ZZ Ω M A ′ ZI Ω M 0 0 0 A ′ IB Ω M A ′ IZ Ω M A ′ II Ω M 0 0 0 0 0 0 A ′ BB Ω CC A ′ BZ Ω CC A ′ BI Ω CC 0 0 0 A ′ ZB Ω CC A ′ ZZ Ω CC A ′ ZI Ω CC 0 0 0 A ′ IB Ω CC A ′ IZ Ω CC A ′ II Ω CC                                   T B Ω M T Z Ω M T I Ω M T B Ω CC T Z Ω CC T I Ω CC                  t+∆t =                  rhs ′ B Ω M rhs ′ Z Ω M rhs ′ I Ω M rhs ′ B Ω CC rhs ′ Z Ω CC rhs ′ I Ω CC                  . (41)
To solve for the interior and boundary values of the mold domain for a time t + ∆t, we have:

A ′ BB Ω M T B Ω M = rhs ′ B Ω M -A ′ BZ Ω M T Z Ω M -A ′ BI Ω M T I Ω M A ′ ZZ Ω M T Z Ω M = rhs ′ Z Ω M -A ′ ZB Ω M T B Ω M -A ′ ZI Ω M T I Ω M . (42)
Since there are 3 unknowns and only 2 equations in Eqn. [START_REF] Storti | PETSc-FEM, A General Purpose, Parallel, Multi-Physics FEM Program, International Center of Computational Method in Engineering (CIMEC)[END_REF], the temperature values at the interpolation boundary are obtained from the interior nodes of Ω CC as follows:

T I Ω M = ΠCC→M T Z Ω CC , (43)
where Π CC→M is an interpolation operator which enables to couple both domains. Then, the full linear system for Ω M is obtained as:

A ′ BB Ω M T B Ω M = rhs ′ B Ω M -A ′ BZ Ω M T Z Ω CC -A ′ ZI Ω M ΠCC→M T Z Ω M , A ′ ZZ Ω M T Z Ω M = rhs ′ Z Ω M -A ′ ZB Ω M T B Ω M -A ′ ZI Ω M ΠCC→M T Z Ω CC . (44)
Analogously, we can obtain a full linear system for Ω CC as follows:

A ′ BB Ω CC T B Ω CC = rhs ′ B Ω CC -A ′ BZ Ω CC T Z Ω CC -A ′ ZI Ω CC ΠM→CC T Z Ω M , A ′ ZZ Ω CC T Z Ω CC = rhs ′ Z Ω CC -A ′ ZB Ω CC T B Ω CC -A ′ ZI Ω CC ΠM→CC T Z Ω M , (45)
where the values at the interpolation boundary T I Ω CC of Ω CC are obtained through the interpolation operator Π M → CC . Then, the following reduced form of the Eqn. ( 41) is obtained:

          A ′ BB Ω M A ′ BZ Ω M 0 A ′ BI Ω M ΠCC→M A ′ ZB Ω M A ′ ZZ Ω M 0 A ′ ZI Ω M ΠCC→M 0 A ′ BI Ω CC ΠM→CC A ′ BB Ω CC A ′ BZ Ω CC 0 A ′ ZI Ω CC ΠM→CC A ′ ZB Ω CC A ′ ZZ Ω CC                     T M B T M Z T CC B T CC Z           t+∆t =           rhs ′ B Ω M rhs ′ Z Ω M rhs ′ B Ω CC rhs ′ Z Ω CC           , (46)
where it can be noticed that both domains are linked by the interpolation operators. Since the resulting systems are not symmetric, their solution can be obtained iteratively with the Bi-Conjugate Gradient Stabilized (BiCGStab) or the Generalized Minimal Residual (GMRes) methods, being the latter the adopted one in this work. The iterative coupling algorithm, the high-order interpolation operator construction, and the coupling accuracy-efficiency are out of the scope of this work, and we refer the reader for a deeper insight to the Refs. [START_REF] Storti | Optimization of an internal blade cooling passage configuration using a chimera approach and parallel computing[END_REF][START_REF] Storti | A matrix-free Chimera approach based on Dirichlet-Dirichlet coupling for domain composition purposes[END_REF][START_REF] Storti | On the performance of a Chimera-FEM implementation to treat moving heat sources and moving boundaries in time-dependent problems[END_REF]. An identical coupling process is performed to transfer information between the Ω P 2 and

the Ω M domains. The coupling between the Ω P and Ω CC domains is not considered in this work since their respective meshes are not expected to face an overlap scenario. In this context, the Ω CC and Ω P domains can be treated as a single mesh in practice, i.e., the domains are meshed independently, and then a single mesh is obtained by renumbering the nodes and connectivities of one of them.

Computational grids for the overlapped domains scheme

Each mesh of the components is conformed by structured tetrahedral linear elements, as shown in Fig.

18

. The Ω P and Ω CC meshes are refined towards the regions of higher thermal gradients, where the ∆x max /∆x min ratios are roughly 34.7 and 28.7, respectively; while the insert mesh is more uniform, being such a ratio of approximately 3.6. A mesh convergence study was performed to assess the proper element sizes needed for the overset scheme to properly capture the temperature field within the part, and to ensure accuracy and efficiency of the coupling. The values of the most suitable multi-domain meshes are depicted in Table 5. It is worth mentioning that the meshes are designed such that the element size ratio at the coupling interfaces is between 1 and 3 to guarantee accuracy of the FEM solution [START_REF] Storti | Optimization of an internal blade cooling passage configuration using a chimera approach and parallel computing[END_REF][START_REF] Storti | A matrix-free Chimera approach based on Dirichlet-Dirichlet coupling for domain composition purposes[END_REF][START_REF] Storti | On the performance of a Chimera-FEM implementation to treat moving heat sources and moving boundaries in time-dependent problems[END_REF].

Domain

Optimization problem description

Mathematical formulation

The mathematical description of the optimization of the CCs' configuration can be defined as:

min x F obj (x) subject to xi L ≤ xi ≤ xi U , for i = 1, ..., n gj(x) = cmin -cj ≤ 0, for j = 1, ..., me, h(x) = T Γ P 1 -T (SI) Γ 1-2 ≤ 0, l(x) = T Ω P 1 -T (SI) Ω 1 ≤ 0, ( 47
)
where x is the vector of the design variables defining the layout of the CC of the fixed and mobile parts of the insert, F obj (x) is the objective function to be minimized, g j (x) is a constraint function which prevents an overset between the CC branches by ensuring a minimum distance c min between its centers c j , h(x) and l(x) are two constraints to ensure that the part and surface temperatures are below the value obtained by the SI (see Section 4.1), and x i L and x i U are the upper and lower bounds of the design variables, respectively (further discussed in Section 5.3.3).

Objective function definition

The desired objective function to minimize comprises two terms, the first one, whose main objective is to minimize the standard deviation of the temperature on the surface of the part, and the second, which aims to reduce the gap between maximum and minimum values of the part temperature which indirectly tends to reduce the hot-spots in the insert. Then, the objective function yields:

F obj (x) = w Γ P 1 T (x, t f ) -T (x, t f ) 2 dΓ F10 + (1 -w) T max Γ P 1 (x, t f ) -T min Γ P 1 (x, t f ) F20 , (48)
where T , T max Γ P 1 and T min Γ P 1 are the average, maximum, and minimum temperatures at the part surface, respectively, F 10 and F 20 are normalization parameters, and w is a weighting factor that can take values in the range [0,1]. In this analysis, the values obtained in case 3 of Table 4 are used for normalization, and a value of w = 0.75 is considered following literature recommendations [START_REF] Agazzi | Contribution à l'optimisation de la thermique des outillages dans le procédé d'injection des polymères thermoplastiques[END_REF][START_REF] Tang | A Computer-Aided Optimization Approach for the Design of Injection Mold Cooling Systems[END_REF].

Design variables of the CCs

The design variable of the optimization problem is conformed by a total of 23 elements. Such variables will define the position of the 16 points shown in Fig. 19, employed for the construction of both CCs' domains, the one that is placed in the fixed part of the insert, and the one placed on its moving part.

Exploiting the fact that the injected part is symmetrical with respect to the yz plane, once the points of only one half of the insert are defined, the others are obtained by performing a proper reflection of the points. Then, the design variables are comprised in the vector x as follows:

x = z 1 ,
where the z i value controls the distance to the cavity of each ith point. Furthermore, the CCs are allowed to bend at the middle (see Fig. 19), and the pitch distances can be modified, which are controlled by the d i parameter defining the position in the x axis of a given CC branch i. Following instructions of the mold manufacturer, in order for the new insert with the optimized CCs be able to fit in the original mold's cavity and mold's core (see Fig. 2), the entrances and exits positions of each CC must remain fixed during the optimization process, thus they are not considered as design variables. The limits of each design variable are defined in Table 6 and were defined according to the next assumptions:

• the z position for each variable considers a minimum distance between the injected part and the CC surface of 9mm, and a maximum distance of 21mm,

• the d position for each variable is thought for the channels to be capable of addressing as much as possible of the insert domain, while avoiding superposition between them, i.e., each branch can move within a restricted region within the insert. The construction of the Ω CC mesh is performed automatically in SALOME [START_REF] Ribes | Salome platform component model for numerical simulation[END_REF] at each OFE with a python script by following Algorithm 1.

ALPSO optimizer

Regarding the optimization task, the ALPSO algorithm [START_REF] Sedlaczek | Using augmented Lagrangian particle swarm optimization for constrained problems in engineering[END_REF] is used in this work, which is able to solve non-linear non-smooth optimization problems that can be subject to equality and/or inequality constraints.

The algorithm is an extended version of the original PSO method developed by Eberhart and Kennedy [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF], with the addition of an augmented Lagrange multiplier approach allowing to deal with constrained problems. [START_REF] Sedlaczek | Using augmented Lagrangian particle swarm optimization for constrained problems in engineering[END_REF]. Specifically, the ALPSO implementation of the open-source software pyOpt [START_REF] Perez | pyOpt: A Python-based object-oriented framework for nonlinear constrained optimization[END_REF] was chosen here due to the following reasons:

• the evolutive-based working principle of PSO makes it more effective than gradient-based optimizers in approaching a global minimum [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF],

• the parallelization capabilities implemented in pyOpt trough the mpi4py [START_REF] Dalcin | Parallel distributed computing using Python[END_REF] library allows to perform a dynamic process management (DPM) to distribute each OFE evaluation on a different core,

• the open license of pyOpt allows its straightforward utilization on a High-Performance Computing (HPC) cluster with the SLURM workload manager [START_REF] Jette | SLURM: Simple Linux Utility for Resource Management[END_REF].

The last point is key due to the amount of OFEs required by stochastic-based optimization algorithms, where (in this work) each one implies large meshes, transient simulations, and several injection cycles.

Particle swarm optimization algorithm

PSO solves problems by considering a population of candidate solutions, namely particles, that move through the search space according to mathematical laws that define their position and velocity. The position and velocity of each particle of the swarm is defined by Eqns. ( 50) and ( 51), respectively, and each constitutes a solution to the problem. First, each particle evaluates the objective function at its random initial position. Then, each particle is updated according to:

p k+1 i = p k i + v k+1 i , (50)
v k+1 i = ωv k i + c 1 r 1 (p p b ,k i -p k i ) + c 2 r 2 (p g b ,k i -p k i ), (51)
where p b stands for the best value that the particle has achieved thus far; g b is the best particle among the entire population; ω is the inertia weight factor; v k i and p k i are the particle velocity and position, respectively, at iteration k; r 1 and r 2 are random numbers between (0,1); and c 1 and c 2 are called cognitive and social parameters.

Optimizer parameters

The convergence performance of stochastic-based optimizers, including the PSO algorithm, is parameterdependent and varies depending on the problem being solved [START_REF] Jiang | Stochastic convergence analysis and parameter selection of the standard particle swarm optimization algorithm[END_REF][START_REF] Shi | Empirical study of particle swarm optimization[END_REF][START_REF] Xueyan | Swarm size and inertia weight selection of particle swarm optimizer in system identification[END_REF]. Based on the study of Xueyan and Zheng [START_REF] Xueyan | Swarm size and inertia weight selection of particle swarm optimizer in system identification[END_REF] a number of 40 particles was selected. Furthermore, the value of ω is set to change from an initial ω 1 to a final ω 2 value during optimization to improve the convergence rate [START_REF] Sedlaczek | Augmented Lagrangian particle swarm optimization in mechanism design[END_REF]. The authors of this work based the adoption of such parameters on the specific literature recommendations and confirmed their suitability through numerical experimentation. The parameters used are depicted in Table 7.

Parameter Description Value

NSW

Overset-based optimization workflow

The optimization workflow is conformed by two blocks (see Fig. 20): a preprocessing stage, which is performed only once previous to the optimization procedure, comprising the CAD generation/manipulation of the insert and the injected part with their corresponding meshing, and an optimization stage which is repeated at each OFE, involving the automatic CAD geometry and meshing generation of the CCs, along with the overset scheme coupled with ALPSO. After the boundary conditions are applied to the B set nodes of each domain, the detection of the interpolation nodes (I set nodes) is performed. Subsequently, the interpolations operators are constructed for the algebraic coupling of the three domains, and the transient solution of the FEM analysis is obtained at each time step for each injection cycle. From the final global solution, the objective and the constraints functions are computed. The values are passed to the optimizer, and, if the stopping criteria is met while accomplishing the constraints, the optimization process is finalized.

Otherwise, a new set of variables is proposed for the CCs' design and the process is repeated. A tolerance value of ϵ = 0.01 and a maximum number of iterations of n max = 600 were chosen as stopping criteria for the present communication.

the mesh would be relegated to the algorithm employed by the meshing software, forcing quality checking before each FEM analysis. In case of changes in the geometry of the part to be injected, it will be required by this approach to remesh solely such component of the computational domain, reducing the meshing efforts. Furthermore, in case a different shape section of the CCs would need to be considered by the mold manufacturer, only the python-based script (see Algorithm 1) to generate the CCs' mesh would need to be 6. Optimization results and discussion

Convergence

The optimization process stopped at iteration 600, reaching the n max value. Even though the convergence tolerance was not reached, the convergence tendency can be clearly observed in Fig. It is observed from Fig. 21 that a smooth and marked convergence rate was obtained by ALPSO coupled with the overset-FEM methodology. As the algorithm started the optimization task by proposing random values of the design variables, ALPSO reached a value of F obj = 1.46. However, as the optimization evolved, a value of F obj = 0.65 was finally achieved, representing a 35% of improvement with respect to the reference case. Furthermore, in Fig. 22 are shown all the trials performed by the optimizer, proving that the stochastic algorithm has widely explored the search space. The tendency to the final optimized configuration can be observed as well.

Optimized CCs' layout analysis

The optimized configuration of the CCs for the fixed and mobile parts of the PI is shown in Fig. 23. It can be seen that the disposition of the optimized CCs is very different from that originally employed in the SI, obtaining a bend configuration where the mid-length of the CCs' branches are closer to the part than the elbows. It is worth noticing that, in order to reduce the thermal gap between maximum and minimum temperatures, the CCs tend to move away from the areas that cool the fastest, as part corners that have large surface-to-volume ratios. The findings of this study align with prior research conducted by Tang et al.

[63] that evaluated the performance of optimized CCs on a 2D ribbed plate, whereby the resulting CCs were situated further away from the corner regions of the plate. Moreover, the results are consistent with those reported by Agazzi et al. [START_REF] Agazzi | Optimal cooling design in injection moulding process -a new approach based on morphological surfaces[END_REF], wherein the optimized solutions indicated that higher coolant temperatures were required in regions of the dilation surface that are proximal to edges, while lower coolant temperatures were necessary for the inner regions of the part.

It can be also observed from Fig. 23 that the pitch distances in the mobile part of the PI remain roughly constant, although such distances are significantly reduced in comparison to the SI CCs' configuration (please see Fig. 23b). Regarding the optimized configuration of the CCs of the insert fixed part, it is observed from Figs. 23a and 23c that the pitch distance is increased at the sharp corner of the lower part.

As aforementioned, due to the higher surface-to-volume ratios, such a region needs to cool down slower than other internal regions in an attempt to achieve homogeneous cooling of the part. It can be further observed from Fig. 23c that two branches of the CCs are close to each other and close to the cavity in the upper region of the injected part. Since the major thickness is located there, it becomes evident from the results that more cooling capacity is required. By comparing the heat flux density at the CCs' surface of the original CCs' disposition (see Fig. 15a) and that of the optimized configuration (Fig. 23), it can be observed that the latter resulted in a qualitatively enhanced heat flow towards the CCs, as compared to the original design.

It was also of interest to compare the optimized solution with that of a non-bent CCs' configuration. to that obtained by the optimized PI, and T max Γ P 1 is approximately 10 • C higher in the plastic mold. Not surprisingly, it is found in the literature that PIs usually lead to higher working temperatures and higher warpage compared to SIs. These results suggest that more research efforts should be put into developing new composite materials capable of reaching higher thermal conductivity values to reduce the cycle time of the IM process while improving the overall thermal performance of the mold.

Convective boundary validation by CFD simulation

This subsection presents the validation process for the heat transfer coefficient estimation obtained from the Dittus-Boelter correlation through the utilization of CFD simulations. The focus of these simulations is to solve for the velocity field inside the pipes using the Navier-Stokes module of PETSc-FEM. Two scenarios are examined: the first scenario encompasses the original arrangement of the CCs in the PI, while the second scenario incorporates the optimized CCs' arrangement.

CFD model

In these 3D simulations, a steady-state condition is assumed for the fluid, which operates under turbulent flow conditions (Re ≈ 16611). To accurately capture the flow behavior, the SST kω turbulence model with the Reynolds-Averaged Navier-Stokes (RANS) approach is employed. It is important to note that the CFD computation is focused exclusively on the reduced computational model, where the velocity and pressure fields are solved solely within the insert domain, as shown in Fig. 25. The SUPG-PSPG (Streamline Upwind Petrov-Galerkin -Pressure Stabilization/Petrov-Galerkin) stabilization technique [START_REF] Garibaldi | Numerical simulations of the flow around a spinning projectile in subsonic regime[END_REF] is used here to cope with the convective-dominant flow problem. Considering that the variations in fluid temperature are not expected to significantly affect the dynamics of the coolant (as elaborated in Section 3.4.3), a weak coupling strategy is implemented. This approach entails solving the velocity and pressure fields independently until convergence is achieved. Subsequently, the obtained velocity field is employed to the stagnation points, is negligible. This is evident from the mean outlet temperatures of 10.016 • C and 10.021 • C for the mobile and fixed CCs, respectively, observed in Fig. 27a. Additionally, it is pertinent to observe that the coolant temperature remains remarkably stable due to the low heat diffusivity exhibited by the insert material, yielding minimal temporal variations. By recalling the heat balance of Eqn. [START_REF] Balthazar | Conception d'un outillage de mise en oeuvre des composites à haute dynamique thermique via l'utilisation de structures lattices[END_REF] and considering t f = 360s, the resulting temperature increase in the coolant is estimated to be approximately ∆T ci-o ≈ 0.032 • C. This estimation is in good agreement with the numerical results obtained. Moreover, the optimized CCs' configuration, characterized by a higher heat flow rate, leads to mean fluid temperatures at the outlet of 10.050 • C and 10.065 • C for the mobile and fixed CCs, respectively. These results align consistently with the estimation provided by Eqn. [START_REF] Balthazar | Conception d'un outillage de mise en oeuvre des composites à haute dynamique thermique via l'utilisation de structures lattices[END_REF], which predicts a temperature increase of ∆T ci-o ≈ 0.049 • C for t f = 240s.

The solutions of the mean temperatures of the part volume, part surface, and cavity surface obtained by the CFD approach are depicted in Fig. 28. The solutions are compared to those obtained by employing the convective boundary condition of Eqn. [START_REF] Feng | Design and fabrication of conformal cooling channels in molds: Review and progress updates[END_REF], with h c = 10100W/(mK) and T c = 10 • C. From the figure, it is evident that both approaches yield negligible temperature differences. Specifically, the mean temperature differences of the cavity surface, part surface, and part volume obtained using both approaches are consistently below 0.2 • C once the PSS is reached in each case. This minimal temperature difference can be attributed to the fact that, as discussed in Section 4.3, the convective heat transfer coefficient does not play a dominant role in the heat transfer from the injected part to the CCs. Instead, the length between the part and the CCs, as well as the thermal conductivity of the PI, are more significant factors. The

Figure 1 :

 1 Figure 1: Detailed scheme showing the part geometry to be injected. Isometric, frontal, and section views are depicted. The A-A section is defined by a cutting plane.

Figure 2 :

 2 Figure 2: Mold of the produced plastic part: (a) injection mold assembly, (b) fixed and mobile parts of the mold, (c) mold's cavity and core with their respective CC and insert halves, and (d) the insert halves with the CC disposition.

Figure 3 :

 3 Figure 3: Cross section of the mold assembly showing the main components and the CCs' disposition.

Figure 5 :

 5 Figure 5: Conformal mesh used for the transient heat transfer FEM analysis.

  b) Conformal approach (this work).

Figure 6 :

 6 Figure 6: Conformal and non-conformal meshes approaches to handle the thermal contact resistances between interfaces.

5 •

 5 C. It can also be observed from Figs. 7c and 7d that the mean temperature of the part is ≈ 208 • C at the end of the cooling stage, while the mean temperature of the cavity surface of the PI is ≈ 175 • C. Evidently, the temperature of the part is very far from the ejection temperature needed for ABS. Furthermore, the insert surface temperature surpasses by 32.8 • C the T g of the PC%10CF (T g = 142.2 • C).

Fig. 9 ,Figure 7 :

 97 Figs.9c and 9b, respectively. For the experimental trials, the cooling time and the coolant temperature were defined in t f = 200s and T c = 20 • C, respectively. As the primary objective of the experimental assessment was to validate the numerical model, rather than optimizing the cycle time to maximize productivity, the ejection times were intentionally prolonged to ensure the attainment of a steady-state condition for each injection cycle. To ensure measurement repeatability, three ABS parts were systematically injected, aiming

Figure 8 :

 8 Figure 8: Temperature field obtained during the injection process of the ABS part: (a)-(b) by employing a steel insert, (c)-(d) by using an identical thermoplastic insert.

Figure 9 :

 9 Figure 9: Multiplate mold: (a) mobile part with the thermoplastic insert detailed in (b), (c) fixed part with the reference steel insert.

Figure 10 :

 10 Figure 10: Comparison of numerical and experimental results. From left to right: an actual photo showcasing the IM process, an IR image capturing the thermal distribution, a postprocessed IR temperature field visualized in Paraview, and the FEM (WE) simulation result also visualized in Paraview. In (a) is depicted the final state of the produced part following the completion of the cooling cycle, while (b) shows the PI cavity immediately after the part ejection. Temperatures are given in • C.

  Temperature at the part surface through two defined lines. Temperature at the cavity surface trough two defined lines.

Figure 11 :

 11 Figure 11: Temperatures comparison between the numerical and the experimental results along two lines at the part (a) and mold's cavity (b) surfaces, at the end of the cooling time.

1 . 5 × 10 - 2 ∼ 1 × 3 PI 6 ∼ 1 . 5 × 10 - 2 ∼ 1 ×

 15102136151021 P-CC /λ 2 R C =1/h c R P =e/(2 1 ) 10 -3 ∼ 5 × 10 -4 ∼ 1 × 10 -4 3.5 ∼ 7.5 × 10 -10 -2 ∼ 6 × 10 -2 ∼ 1 × 10 -4 3.5 ∼ 7.5 × 10 -3 (b)Table of order-of-magnitude.

Figure 12 :

 12 Figure 12: Overall stationary heat transfer electrical analogy between the core of the injected part and the CC surface. e is in mm and the thermal resistances in m 2 K W -1 .

  0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 conductivity [W/(mK)] Cycle time [s] PI: LP-CC = 21mm, e = 3.50mm PI: LP-CC = 21mm, e = 4.75mm PI: LP-CC = 21mm, e = 6.00mm PI: LP-CC = 9mm, e = 3.50mm PI: LP-CC = 9mm, e = 4.75mm PI: LP-CC = 9mm, e = 6.00mm SI: LP-CC = 21mm, e = 3.50mm SI: LP-CC = 21mm, e = 4.75mm SI: LP-CC = 21mm, e = 6.00mm

Figure 13 :

 13 Figure 13: Comparison of injection process cycle times for varying PI thermal resistance values and three different thicknesses (e) of the injected part, in comparison to the cycle times achieved with a SI. SI cycle times correspond to fixed material

  resistance [m 2 K W -1 ] Cycle time [s] -Plastic insert Thermal contact resistance [m 2 K W -1 ] Cycle time [s] -Steel insert PI: e = 3.5mm PI: e = 4.75mm PI: e = 6mm SI: e = 3.5mm SI: e = 4.75mm SI: e = 6mm

Figure 14 :

 14 Figure 14: Cycle time required to achieve the ejection temperature of the part as a function of the thermal contact resistance.

Figure 15 :

 15 Figure 15: Results of temperature distribution at the part surface and heat flux density on the CC surfaces at the end of the cooling stage for three positions of the CC.

Figure 16 :

 16 Figure 16: Temperature distribution on the part surface at the end of the cooling stage.

  Ω CC Interpolation nodes (I) of Ω M Interpolation nodes (I) of Ω P Boundary nodes (B) of Ω M Boundary nodes (B) of Ω CC (b) Domains structured meshing and nodes splitting.

Figure 17 :

 17 Figure 17: Domain splitting framework for the overset Chimera scheme. In the detailed view the nodes splitting of each subdomain. The remaining nodes are those of the set Z (interior nodes) and are not shown here.

Figure 18 :

 18 Figure 18: Structured meshes of the Ω P , Ω M , and Ω CC domains used for the optimization procedure with the overset scheme. At the left of the figure, the three domains overlapped, and at the middle the Ω P (top) and Ω CC (bottom) meshes with their respective detailed views at right. The background Ω M mesh is practically uniform, and the Ω P and Ω CC are anisotropically refined towards the regions of highest thermal gradients.

Figure 19 :

 19 Figure 19: Variables employed to construct the CCs' domain. The variables are designed to always consider a minimum of 9mm distance between the part-cavity and the CC surface.

Algorithm 1 :

 1 Automatic construction of a structured 3D mesh for cooling channels Data: Vertex positions given by the optimizer Result: Structured 3D-mesh of cooling channels to be used by overset-FEM Input: Vector with optimizer variables x Output: xnod,icone,I,B,Z 1 use x to define vertex positions in the 3D space 2 reflect points over yz plane 3 construct wire to define the axis of the CCs 4 perform fillet operations on the wire to smooth the CC curves 5 perform structured 2D meshes at the entrances 6 sweep the 2D mesh through wire obtained in point 4 → hexahedral structured 3D mesh 7 split hexahedral elements into tetrahedrons 8 define node groups B, I, Z 9 export: nodes coordinates (xnod), elements connectivities (icone), and groups B, I and Z

Figure 21 :

 21 Figure 21: Convergence of ALPSO with the overset methodology. The values obtained by the PI case 3 with L P -CC = 9mm are used for normalization.

Figure 22 :

 22 Figure 22: Configurations of the CCs proposed by ALPSO (black thin lines) and the best configuration achieved (red thick lines).

Consequently, incorporating an50Figure 24 :

 24 Figure 24: Temperature distribution on the part surface at the end of the cooling stage obtained by the optimized CCs' configuration.

Figure 25 :

 25 Figure 25: Boundary conditions and structured mesh employed in the CFD validation.

  (a) Streamlines of the velocity with the original CCs. (b) Streamlines of the velocity with the optimized CCs.

Figure 26 :

 26 Figure 26: CFD analysis of the cooling fluid inside the channels of the fixed and mobile parts of the insert, for both configurations, the reference and the optimized one.

Figure 28 :

 28 Figure 28: Comparison of the results obtained by using a convective heat transfer coefficient (hc) with respect to the results obtained by using CFD analysis of the cooling fluid inside the channels. The original (a) and optimized (b) configurations were assessed.

Figure B. 1 :

 1 Figure B.1: Time step adoption.

Figure C. 1 :hc- 1 )ΩF10

 11 Figure C.1: Reflectivity values obtained by the Bruker 80V device for the ABS and PC10%CF samples. The emissivity values were obtained from the assumption of the samples as opaque bodies, and by applying the Kirchoff's law for thermal radiation.

Table 1 :

 1 Material properties of the multi-plate mold and the injected polymer.

	1030	1800	1.07 × 10 -7

Table 3 :

 3 Comparisons

		36.68	30.22	3.35	52.03	39.73	7.17
	FEM (WE)	38.95	32.30	4.23	55.71	42.51	9.71
	FEM (WOE)	39.26	32.96	4.17	56.07	45.36	8.16
	Diff. IR-FEM (WE)	2.27	2.08	0.88	3.68	2.78	2.54
	Diff. IR-FEM (WOE)	2.58	2.74	0.82	4.04	5.63	0.99

of the results obtained by the half-steel and half-plastic insert mold by the experiments compared to numeric results.

Table 4 :

 4 Results obtained by the parametric study.

Table 5 :

 5 Details of the meshes employed for the overlapped scheme.

  z 2 , ..., z 6 , d 1 , d 2 , z 8 , ..., z 16 , d 3 , d 4 , ..., d 7

	Mobile part CCs	, z 7 Fixed part CCs	,

Table 6 :

 6 Variables range for the cooling channels optimization of the fixed and mobile parts of the insert. The distances are in mm.

Table 7 :

 7 ALPSO optimizer parameters employed in pyOpt.

  [START_REF] Reyes | Inverse heat transfer optimization of stamping with overmolding process involving high performance thermoplastic composites: experimental validation[END_REF]. Each OFE took approximately 8hs of computation time on one single core, where on average 6-8 injection cycles must be performed to reach the PSS in each case. The values of ∆t f and ∆t e were set in 5s. The optimization process was parallelized on 40 cores from the computer cluster of Le Centre de calcul intensif des Pays de la Loire, and the entire optimization solution took roughly 5 days.

		1.60
					PI case 1
		1.40			PI case 2 PI case 3 (normalization)
					ALPSO trials
		1.20			Optimized case
	F obj (x)	1.00
		0.80
		0.60
		0	100	200	300	400	500	600
			Objective function evaluations

assumption between boundaries, the heat flux released by the injected part through its surface is modeledΩ 1 Ω 2 Ω 3 Ω 4 Γ 4(out) Γ 4(cc) Γ 3(cc) Γ 2(cc) Γ 3(out) Γ 3-4 Γ 2-4 Γ 1-2 Γ 2-3Figure4: Two-dimensional scheme of the computational domains with their respective boundaries considered for thermal analysis.

i=1 Γ2-i 1 T CR 2-i N T NdΓ,(21)

i=2 Γ3-i 1 T CR 3-i N T NdΓ,(22)

i=2 Γ4-i 1 T CR 4-i N T NdΓ.(23)Finally, the following load vectors related to the heat flux transferred to the CCs/ambient are obtained:F Ω2 = -Γ 2(CC) h c T c N T dΓ,(24)F Ω3 = -Γ 3(CC) h c T c N T dΓ -Γ 3(ext) h eq T env N T dΓ,(25)

Acknowledgments

The authors gratefully acknowledge the financial support received from the Banque Publique d'Investissement (BPI) for the project SAMFAST: Fabrication additive rapide d'empreintes d'outillages d'injection plastique.

The support given in the framework of the project by the enterprises CERO, NANOVIA, IRMA and COM-POSITIC is also gratefully acknowledged.

Coupling of the heat transfer algebraic systems

Once the nodes of each domain are split into groups, a coupling strategy must be followed for transferring information between domains to achieve a final solution of the system. To do so, consider Eqn. [START_REF] Perez | pyOpt: A Python-based object-oriented framework for nonlinear constrained optimization[END_REF] only The fact that the Ω P and Ω M domains do not change throughout the optimization procedure proposed in this work allows for the construction of high-fidelity structured meshes that are maintained on each OFE.

Otherwise, a complete remeshing would be required on each OFE by conventional FEM approaches. Besides the additional computational cost implied in remeshing the entire domain for fine meshes, the quality of changed.

order of several millimeters can significantly impact the mold's thermal performance. It is noteworthy that the enhancement in heat transfer achieved by the proposed scheme may be constrained by the part's geometry, suggesting the potential for further improvements in parts with intricate surfaces or deeper concavities. Further quantitative results are presented in Table 8 where it is observed that for the same IM time parameters of the PI (case 3), i.e. t f = 240s and t c = 300s, the optimized configuration was able to obtain a reduction of the T Γ P 1 of almost 5 solve the transient convection-diffusion equations in the fluid (Ω F ), insert (Ω M ), injected part (Ω P ), and cooling channel domain (Ω CC ). A zero velocity condition is imposed in the latter three domains.

Fluid domain discretization and boundary conditions

A structured mesh was utilized to discretize the fluid domain, comprising a total of approximately 6.240.514 tetrahedral elements. Specifically, Ω F of the mobile part of the mold was represented by 2.821.084 elements, while Ω F of the fixed part was discretized with 3.419.430 elements. The mesh configuration aimed to ensure precise computation of the velocity field and pressure losses, with careful consideration given to positioning the nodes closest to the wall within the viscous sublayer (y + < 5) to obtain accurate solutions.

Boundary conditions applied to CCs included a flow rate of V F = 9lmin -1 at an inlet temperature of

At the surface of the CCs, a non-slip boundary condition was enforced. The outlet was characterized by a zero pressure condition and an adiabatic condition. The remaining thermal boundary and initial conditions corresponded to the specific PI cases, designated as PI: case 1 and PI: case opt., as outlined in the previous sections. Time steps of ∆t f = 5s and ∆t e = 5s were employed. The physical properties of the coolant adopted in the simulations comprised ρ F = 999.75kg m -3 , Cp F = 4200J/(kgK), µ F = 1.30 × 10 -3 Pa s, and κ F = 0.578W/(mK), representing the density, heat capacity, dynamic viscosity, and thermal conductivity of water, respectively.

CFD results

The streamlines of the velocity field derived from the CFD simulations are depicted in Fig. 26

Concluding remarks and future work

This work began with precise 3D finite element transient simulations to assess the feasibility of replacing a traditional steel insert, featuring conformal cooling channels, with a thermoplastic alternative manufactured via fused deposition modeling. Experimental validation through thermal imaging of a half-steel, half-plastic mold scenario allowed us to validate the numerical model, achieving minimal temperature discrepancies. Simulations results under the steel-case processing parameters showed that, despite the usage of cooling channels, the mold's cavity surpassed by 32.8 • C the glass transition temperature of the carbon fiber-reinforced polycarbonate, risking mold failure and preventing ABS polymer ejection. Consequently, modifications to the process parameters and the cooling channel design framework were mandatory, resulting in the following main findings from parametric studies on the plastic mold:

• High sensitivity of cycle times to mold thermal conductivity was observed in the 0.1 -1.0W/(mK) range, e.g., a 0.1W/(mK) polymer resulted in a cycle time 22.7 times longer compared to steel, while the 0.33W/(mK) composite reduced it to 6.4 times.

• Polymer-cavity thermal contact resistances below 10 -3 m 2 K/(W) have a negligible impact on the cycle time for the polymeric insert, while this threshold shifts to 10 -4 m 2 K/(W) for the steel counterpart.

• Decreasing the cavity-to-CCs distance negatively impacted mold thermal performance, increasing the maximum temperature at the part's surface and potentially leading to hot spots.

• Reduced heat dissipation led to smaller through-thickness thermal gradients, with 2-5°C differences in surface and volume temperatures, compared to 21°C with the steel insert.

• Regardless of cavity-to-CCs distance, increased non-uniformity at the part's surface was observed, with standard deviations of roughly 11 • C, in contrast to 2.3 • C obtained with the steel reference insert.

Finally, the proposed scheme to optimize the CCs' layout based on overset meshes coupled with the particle swarm-based optimizer showed favorable convergence rates, resulting in a 42% reduction in the temperature deviation term, a decrease of more than 6 • C in the maximum temperature difference on the part's surface, and a nearly 5 • C reduction in its mean temperature. Accurate CFD simulations validated the bent CCs' configuration obtained by the proposed optimization framework. Future work will focus on manufacturing and experimentally studying the thermal performance of the plastic insert with the optimized CC configuration, as well as on enhancing the thermal conductivity of the insert material.

Appendix A. Mesh convergence and analytical comparison

Appendix A.1. Mesh convergence For mesh convergence analysis, the cyclic 3D transient heat transfer problem was solved on three meshes (M1;M2;M3) with different refinement levels. Anysotropically refined unstructured meshes with tetrahedral elements were constructed with SALOME [START_REF] Ribes | Salome platform component model for numerical simulation[END_REF]. The meshes were constructed so as to have a refinement towards the regions with higher gradients, such as the part-mold, mold-mold, and CCs-mold interfaces. In Table A.1, the number of cells of each mesh, the maximum ∆x max and minimum ∆x min element sizes, the mean temperature of the part surface T Γ1-2 and volume T Ω1 , and the terms accounting for temperature uniformity T dev Γ1-2 and temperature gap ∆T at the part surface are reported. In addition to the data presented in Table A.1, the mean temperatures of the insert surface T Γ2-1 and volume T Ω2 , and part surface and volume T Ω1 , are plotted in Figs. A.1a and A.1b. It is observed from the results that the M1 mesh overestimates the insert temperatures, while the part temperatures are underestimated for such a mesh. The results improve as the mesh is refined. Finally, and for the sake of accuracy, the M3 mesh is adopted for the analyses.

where ρ p , Cp p and e are the density, heat capacity and thickness of the injected part, κ m denotes the thermal conductivity of the mold, and L D and L P -CC are the pitch and part-to-CCs distances, respectively.

By considering the thermophysical properties of the ABS injected polymer and the steel insert reported in Table 1, the IM process parameters of Table 2, the geometric parameters L P -CC = 21mm, L D = 40.5mm, d = 9mm, and T eje = 95.7 • C, then the following values are obtained for three thicknesses of the part:

• for e = 3.5mm : T m = 69.13[ • C],

• for e = 4.75mm :

• for e = 6mm : T m = 72.08[

Appendix C. Emissivity measurements

The radiative properties of the composite (PC10%CF) and ABS samples were measured using a Bruker 80V Fourier Transform spectrometer, in accordance with standard laboratory protocols. This device is capable of combining different sets of sources, separators, and detectors to cover the spectral range of interest. Normal-hemispherical reflectance measurements were performed using integrating spheres, and the absorbance and normal emissivity were deduced indirectly using the first and second laws of Kirchhoff for radiation. The measurements were performed at ambient temperature. Nitrogen was used for cooling the IR detectors. The measurements were carried out on a printed composite PC10%CF sample, and one ABS sample obtained from a manufactured part. A 2mm beam diameter was employed for the measurements.

The results for a range of wavelength of 4µm -16µm are depicted in