<figure><img src="image_1.png" /></figure>

HAL
open science

# On the use of reduced bases in optimization of injection molding

Francisco Chinesta, Fabrice Schmidt

► To cite this version:

Francisco Chinesta, Fabrice Schmidt. On the use of reduced bases in optimization of injection molding. José Antonio Covas; Antonio Gaspar-Cunha. Optimization in polymer processing, Nova Science Publishers, pp.169-194, 2011, 9781611228182. hal-01007527

HAL Id: hal-01007527

https://hal.science/hal-01007527

Submitted on 11 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the use of reduced bases in optimization of injection molding

**Francisco Chinesta**

EADS Corporate Foundation International Chair

GEM UMR CNRS – Centrale Nantes, France

**Fabrice Schmidt**

Institut Clément Ader

ENSTIMAC, Ecole de Mines d'Albi, France

## 1. Introduction

### 1.1. Process description – Injection molding

About 30% of the annual polymer production is transformed by injection molding. It is a cyclic process of forming a plastic into a desired shape by forcing the molten polymer under pressure into a hollow cavity (Agassant et al, 1991; Osswald, 1998). For thermoplastics polymers, the solidification is achieved by cooling. Typical cycle times range from 1 to 100 seconds and depend mainly on the cooling time. The complexity of molded parts is virtually unlimited, sizes may range from very small (<1mm) to very large (>1m), with an excellent control of tolerances.

### 1.2. The injection molding equipment

The reciprocating screw injection molding machine is the most common injection unit used (figure 1). These machines consist of two basic parts, an injection unit and a clamping unit. The injection unit melts the polymer resin and injects the polymer melt into the mold. Its screw rotates and axially reciprocates to melt, mix, and pump the polymer. A hydraulic system controls the axial reciprocation of the screw, allowing it to act like a plunger, moving the melt forward for injection. The clamping unit holds the mold together, opens and closes it automatically, and ejects the finished part.
<figure><img src="image_3.png" /><figcaption>**Figure 1** - Injection molding machine with reciprocating screw (Wesselmann, 1998)</figcaption></figure>

**1.3. Description of the injection molding cycle**

The process is started by plasticizing the material. A resin supplied in the form of pellets or powder is fed from the hopper into the injection unit which consists of a reciprocating screw in a barrel heated by several heating bands. The rotating screw transfers the solid material towards the heated zones of the barrel. The granules melt under the combined action of the heater bands and the dissipated energy induced by shearing of the viscous material when rotating the screw.

The screw stops rotating when the required amount of material has been dosed. In the soak time until injection, the polymer melts by heat conduction from the barrel to the polymer. Before injection we ideally obtain a completely molten, low viscous material with a homogenous temperature distribution.

The injection takes only a small portion of the cycle time. The screw acts like a plunger and pushes the melt into the mold. The polymer flows from the nozzle to the mold which is coupled to the nozzle by a sprue bushing. The melt flows to a cavity by runners and is fed to the cavity through a gate. The gate is simply a restriction in the flow path just ahead of the mold cavity and serves to direct the flow of the melt into the cavity and to limit back flow.

The cooling of the melt starts when the melt leaves the heated section and gets in contact with the cooled mold walls. During filling the hydraulic pressure on the screw is adjusted to follow a programmed transitional speed profile which allows controlling the flow front speed in the cavity.
When the cavity is completely filled the polymer pressure increases instantly. The machine must now instantly stop pushing forward the screw to avoid an over packing of the cavity.

The specific volume of thermoplastic polymers reduces when passing from the molten to the solid state. More melt must be added to the cavity during solidification to compensate for the contraction of the polymer. Therefore, a constant hydraulic pressure is applied so that the screw holds the melt under pressure and pushes more melt into the mold. Eventually, the plastic in the gate freezes, isolating the mould from the injection unit. While the part continues to cool, the melt for the next shot is dosed. At the end of the cooling time the part must be solid enough to retain the shape given by the cavity. The clamping unit opens the mould, ejects the part and closes and clamps the mould again. The next injection cycle starts.

## 1.4. Importance of cooling step for manufacturing injected parts

Part cooling during injection molding is the critical step as it is the most time consuming. An inefficient mold cooling may have dramatic consequences on cycle time and part quality and may require expensive mold rectification. Depending on the wall thickness of the molded parts, it usually takes the major portion of the cycle time to evacuate the heat. Polymers are bad conductors, as thermal conductivity ranges form 0.1 W/mK to 1.8 W/mK (Osswald, 1998). The cooling cycle can represent more than 70% of the injection cycle (Opolski et al, 1987; Lu, 1996). The cooling rate is an important factor for productivity, and important benefits can be achieved by decreasing the cooling time of parts with hot zones badly cooled. A bad design of the cooling channels may generate zones with higher temperatures in the mold, increasing the cooling time.

In addition, different types of injection defects due to a bad thermal regulation of the mold can appear: dimensional defects, structural defects and aspect defects (Yang et al, 1996; Moller et al, 1998).

In order to reduce mold and production costs, an automatic optimization of the geometry of the cooling device and processing parameters (temperature, flow rate...) may be developed. The optimization procedure necessitates to compute numerous transient heat balance problems (eventually non linear). For solving the thermal problem, we need an efficient meshing technique. The Boundary Element Method (BEM) is well adapted for such a problem, because it only requires a surfacing mesh. The displacement of the cooling channels after each optimization iteration is then facilitated (no remeshing). In addition, reduced modeling is useful in order to reduce the CPU time of the direct computations, particularly for 3D computations. Injection molding is a cycle process; which imply the computation of numerous cycles. On figure 2, an example of temperature history of 40 cycles is plotted versus time. Curves (a) and (b) give respectively the maximum of the temperature in the cavity before and after optimization. Curve (c) represents the average temperature at the cavity surface after optimization.
<figure><img src="image_4.png" /><figcaption>**Figure 2** -Temperature history of the first 40 cycles.</figcaption></figure>

# 2. Mold cooling optimization

## 2.1. Introduction

Several CAD and simulation tools are available to help designing the cooling system of an injection mold. Simulation of heat transfer during injection can be used to check a mold design or study the effect of a parameter (geometry, materials...) on the cooling performance of the mold. Several numerical methods such as Finite Element Method (FEM) (Boillat et al, 2002) or Boundary Element Method (BEM) (Bialecki et al, 2002) can be used. Bikas et al (1999) used C-Mold® simulations and design of experiments to find expressions of mean temperature and temperature variation as functions of geometry parameters of the mold.

Numerical simulation can also be used to perform an automatic optimization of mold cooling. Numerical simulation is used to solve the thermal equations and evaluate a cost function related to productivity or part quality. An optimization method is used to modify the parameters and improve the thermal performance of the mold. Tang et al (1998) used 2D transient FEM simulations coupled with Powell's optimization method (Fletcher, 1987) to optimize the cooling channel geometry to get uniform temperature in the polymer part. Huang et al (2001) used 2D transient FEM simulations to optimize the use of mold materials according to part temperature uniformity or cycle time. Park et al (1998) developed 2D and 3D stationary BEM simulations in the injection molds coupled with 1D transient analytical computation in the polymer part (throughout the thickness). The heat transfer integral equation is differentiated to get sensitivities of a cost function to the parameters. The calculated sensitivities are then used to optimize the position of linear cooling channels for simple shapes (sheet, box).

In the next section, we present the use of the Boundary Element Method (BEM) and DRM applied to transient heat transfer of injection moulds. The BEM software, (Mathey et al, 2004; Pirc et al, 2009), was combined with an adaptive reduced
modeling (described extensively later). This procedure will be fully described in section 3, it allows reducing considerably the computing time during the linear system solution in the transient problem. Then, we present a practical methodology to optimize both the position and the shape of the cooling channels in injection molding processes (section 4). We couple the direct computation with an optimization algorithm such as SQP (Sequential Quadratic Programming).

## 2.2. BEM for transient heat balance equation

Using BEM, only the boundary of the domain has to be meshed and internal points are explicitly excluded from the solution procedure. An interesting side effect is the considerable reduction in size of the linear system to be solved (Pirc et al, 2009). The transient heat conduction in a homogeneous isotropic body $\Omega$ is described by the diffusion equation, where $\alpha$ is the material thermal diffusivity, assumed to be constant:

$$ \Delta T(\mathbf{x}, t) = \frac{1}{\alpha} \frac{\partial T(\mathbf{x}, t)}{\partial t} \qquad (1) $$

We define the initial conditions and the boundary conditions (figure 3) as:

$$ \left\{ \begin{array}{l} T(\mathbf{x}, t=0) = T^o(\mathbf{x}) \\ -\lambda \nabla T \cdot \mathbf{n} = q_p \quad \forall \mathbf{x} \in \Gamma_p \\ -\lambda \nabla T \cdot \mathbf{n} = h_c (T - T_c) \quad \forall \mathbf{x} \in \Gamma_c \\ -\lambda \nabla T \cdot \mathbf{n} = h_a (T - T_a) \quad \forall \mathbf{x} \in \Gamma_M \end{array} \right. \qquad (2) $$

Where $\lambda$ is the thermal conductivity (the medium is assumed as homogeneous and isotropic), $\Gamma_p$ is the boundary of the cavity surface (plastic part), $\Gamma_c$ the boundary of the cooling channels and $\Gamma_M$ the mold exterior surface. The temperature of the coolant is $T_c$ and the heat transfer coefficient, $h_c$, represents the heat transfer coefficient between the mold and the coolant. $h_a$ represents the heat transfer coefficient between the mold and the ambient air at a temperature $T_a$. In order to avoid multi-domains calculation and save computation time, the plastic part is taken into account via a heat flux $q_p$ imposed on the mold cavity surface. The flux density $q_p$ is calculated from the cycle time and polymer properties (Mathey et al, 2004; Pirc et al, 2008).
<figure><img src="image_5.png" /><figcaption>Figure 3 –Boundary conditions applied to the mold</figcaption></figure>

Different strategies are possible to solve such problems using BEM. Pasquetti et al (1995) propose to use space and time Green's function. Another solution should be to apply Laplace (Sutradhar et al, 2002) or Fourier (Godinho et al, 2004) transforms on time variable before spatial integration. In this section, we will use only space Green's function inspired on the stationary heat transfer problem (i.e. Laplace's equation)

The basic steps are in fact quite similar to those used for the finite element method. We firstly must form an integral equation from the previous equation by using a weighting integral equation and then use the Green-Gauss theorem:

$$ \int_{\Gamma} T^{*} q d\Gamma - \int_{\Omega} \nabla T^{*} \nabla T d\Omega = \int_{\Omega} \frac{1}{\alpha} \frac{\partial T}{\partial t} T^{*} d\Omega \quad (3) $$

where $T^*$ is the weighting function and $q$ the normal temperature gradient. This is the starting point for the finite element method. To derive the starting equation for the boundary element method, we use the Green-Gauss theorem again to the second left-hand integral. This gives:

$$ \int_{\Gamma} T^{*} q d\Gamma - \int_{\Gamma} T q d\Gamma + \int_{\Omega} \Delta T^{*} T d\Omega = \int_{\Omega} \frac{1}{\alpha} \frac{\partial T}{\partial t} T^{*} d\Omega \quad (4) $$

For the boundary element method we choose $T^*$ to be the fundamental solution of Laplace's equation (also called the Green' function):

$$ \Delta T^* + \delta_x = 0 \quad (5) $$

where $\delta_x$ is the Dirac delta distribution at the point $x$ located inside the domain $\Omega$ or on its surface $\Gamma$. If we use now the integral property of the Delta function, we obtain (Brebbia et al, 1992):
$$ \int_{\Omega} \delta_x T \, d\Omega = c(\mathbf{x})T(\mathbf{x}) \qquad (6) $$

with:

$$ c(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega \\ \frac{1}{2} & \text{if } \mathbf{x} \in \Gamma \text{ and } \Gamma \text{ is smooth at } \mathbf{x} \\ \frac{\text{internal angle}}{2\pi} & \text{in 2D \& if } \mathbf{x} \in \Gamma \text{ (}\Gamma \text{ is not smooth at } \mathbf{x}) \\ \frac{\text{inner solid angle}}{2\pi} & \text{in 3D \& if } \mathbf{x} \in \Gamma \text{ (}\Gamma \text{ is not smooth at } \mathbf{x}) \end{cases} \qquad (7) $$

Thus, we get the integral equation:

$$ c(\mathbf{x})T(\mathbf{x}) + \int_{\Gamma} Tq^* d\Gamma + \int_{\Omega} \frac{1}{\alpha} \frac{\partial T}{\partial t} T^* d\Omega = \int_{\Gamma} T^* q d\Gamma \qquad (8) $$

The fundamental solutions of Laplace's equation are well-known (Brebbia et al, 1992). $T^*$ and $q^*$ are then defined by equations (9) and (10) depending on the dimension of the problem:

$$ \text{For 2D problems: } \begin{cases} T^* = -\frac{1}{2\pi} \ln(r) \\ \frac{\partial T^*}{\partial n} = \frac{d}{2\pi r^2} \end{cases} \qquad (9) $$

$$ \text{For 3D problems: } \begin{cases} T^* = \frac{1}{4\pi r} \\ \frac{\partial T^*}{\partial n} = \frac{d}{4\pi r^3} \end{cases} \qquad (10) $$

where $r = \|\overline{PM}\|$ and $d = -\vec{r} \cdot \vec{n} = -\mathbf{r} \cdot \mathbf{n}$ (figure 4):
<figure><img src="image_6.png" /><figcaption>**Figure 4** –Definition of distances used to compute fundamental solutions</figcaption></figure>

Again, a similar step to FEM consists in meshing. So, the boundary Γ of the domain is divided into N_c elements.

To express now the domain integral in terms of equivalent boundary integrals, we introduce the DRM approximation (Brebbia et al., 1999). The DRM consists in seeking the solution as a series of particular solutions $\hat{T}$ and $\hat{q}$ interpolated on N points inside and on the boundary of the domain.

$$ \Delta T = \sum_{k=1}^{N} \alpha_k \Delta \hat{T}_k = \text{with } \alpha = \frac{1}{\alpha} \mathbf{F}^{-1} \mathbf{D} \qquad (11) $$

where the vector $\alpha$ contains the coefficients $\alpha_k$, vector $\mathbf{D}$ the different derivatives $\frac{\partial T_k}{\partial t}$ and matrix $\mathbf{F}$ consists of the values of the interpolation function at each point.

A commonly used interpolation function is the polynomial radial function leading to particular solutions (Brebbia et al., 1992). Applying BEM to the modified equation (11) leads to the new linear system of equations (12).

$$ \mathbf{H}\mathbf{T} - \mathbf{G}\mathbf{Q} = \frac{1}{\alpha}(\mathbf{H}\hat{\mathbf{T}} - \mathbf{G}\hat{\mathbf{Q}})\mathbf{F}^{-1}\mathbf{D} \qquad (12) $$

where $\mathbf{H}_{ij} = c\delta_{ij} + \int_{\Gamma_i} q^* d\Gamma$ and $\mathbf{G}_{ij} = \int_{\Gamma_i} T^* d\Gamma$. A Newmark time scheme is applied to the temperature and flux leading to equation (13):

$$ \left( \theta \mathbf{H} + \frac{1}{\Delta t} \mathbf{C} \right) \mathbf{T}^{p+1} - \theta \mathbf{G} \mathbf{Q}^{p+1} = - \left( (1-\theta) \mathbf{H} - \frac{1}{\Delta t} \mathbf{C} \right) \mathbf{T}^p + (1-\theta) \mathbf{G} \mathbf{Q}^p \qquad (13) $$
where $C = -\frac{1}{\alpha}(\hat{H}\hat{T} - \hat{G}\hat{Q})$, $\Delta t$ is the time step and from the chosen values of $\theta \in [0,1]$ the time integrating scheme will result.

## 2.3. Coupling DRBEM with an optimization method

The heat transfer computation is coupled with an optimization method to automatically modify the parameters at each optimization iteration as shown in figure 5. Given the initial parameters, the DRBEM simulation is performed and the cost function is calculated. The optimization method allows updating parameters according to constraints until a minimum of the cost function is found. SQP (Sequential Quadratic Programming) (Mathey, 2004) is used for the optimization of continuous non linear functions with continuous non linear constraints.

<figure><img src="image_7.png" /><figcaption>**Figure 5** – Heat transfer simulation / optimization coupling</figcaption></figure>

In order to reduce the computing time during the linear system solution (especially in 3D), we propose the use of the model reduction within the BEM solver. Before presenting the results of optimization, the reduced modeling approach is summarized.

# 3. Reducing modeling

## 3.1. Introduction

Most engineering systems can be modeled by a continuous model usually expressed by a system of linear or non-linear coupled partial differential equations describing
the different conservation balances (momentum, energy, mass and chemically reacting substances). From a practical point of view, the determination of its exact solution, that is, the exact knowledge of the different fields characterizing the physical system at any point and time instant (velocity, pressure, temperature, chemical concentrations, ...), is not possible in real systems due to the complexity of models, geometries and/or boundary conditions. For this reason, the solution is searched only at some points and at some times, from which it could be interpolated to any other point and time. Numerical strategies allowing this kind of representation are known as discretization techniques. There exist numerous discretization techniques, e.g. finite elements, finite volumes, boundary elements, finite differences, meshless techniques, among many others. The optimal technique to be applied depends on the model and on the domain geometry. Progresses in numerical analysis and in computation performances make currently possible the solution of complex systems involving millions of unknowns related to the discrete model. However, the complexity of the models is also increasing exponentially, and today engineers are not only interested in solving a model, but also in solving these models many times (e.g. when they address optimization or inverse identification). For this purpose, strategies able to speed-up the numerical solution, preserving the solution accuracy, are in focus.

In the context of control, optimization or inverse analysis, numerous problems must be solved, and for this reason the question related to the computation time becomes crucial. The question is very simple: is it possible to perform very fast and accurate simulations? Different answers have been given to this question depending on the scientific community to which this question is addressed. For specialists in computational science the answer to this question concerns the improvement of computational resources, high performance computing and the use of parallel computing platforms. To some specialists in numerical analysis, the challenge is in the fast resolution of linear systems via the use of preconditioners or multigrid techniques, among many others. To others, the idea is to adapt the cloud of nodes (points where the solution is computed) in order to avoid excessive number of
unknowns. Many other answers have been given, however at present all these approaches allow to slightly alleviate the computation efforts. However, the fast and accurate computation remains a real challenge.

This section describes a different approach based on model reduction, allowing fast and accurate computations. The idea is very simple. Consider a domain where a certain model is defined, as well as the associated cloud of nodes able to represent by interpolation the solution everywhere. In general, the number of unknowns scales with the number of nodes and, for this reason even if the solution is evolving in time smoothly all the nodes are used to describe it at each time step. In the reduced modeling that we describe later the numerical algorithm is able to extract the optimal information describing the evolution of the solution in the whole time interval. Thus, the evolution of the solution can be expressed as a linear combination of a reduced number of functions (defining the reduced approximation basis), and then the size of the resulting linear problems is very small, and consequently the CPU time savings can attain several orders magnitude (sometimes, in the order of millions).

The extraction of this relevant information is a well known topic based on the application of the proper orthogonal decomposition, also known as Karhunen-Loève decomposition (Karhunen, 1946; Loève, 1963) that is summarized in the next section. This kind of approach has been widely used for weather forecast purposes (Lorenz, 1956), turbulence (Sirovich, 1987; Holmes et al., 1997), solid mechanics (Krysl et al., 2001) but also in the context of chemical engineering for control purposes (Park and Cho, 1996).

Usually reduced modeling performs the simulation of a similar problem or of the desired one in a short time interval. From these solutions the Karhunen-Loève decomposition applies, allowing the extraction of the most relevant functions describing the solution evolution. Now, it is assumed that the solution of a “similar” problem can be expressed using this reduced approximation basis, allowing for a significant reduction on the discrete problem size and then to significant CPU time
savings. However, in general, the question related to the accuracy of the computed solutions is ignored. An original approach combining the model reduction and the control of the solution accuracy was proposed by Ryckelynck (2005), and applied later in a large catalogue of applications (Ryckelynck et al., 2006; Ryckelynck et al., 2005; Ammar et al., 2006, Niromandi et al., 2008; Chinesta et al., 2008; Verdon et al., 2009). This model reduction strategy can be coupled with usual finite element or boundary element discretizations (Ammar et al., 2009).

We summarize in this section the main ideas of this reduction strategy for the non specialist in numerical analysis, in order to show its potentiality in many domains of engineering and, in particular, in the context of optimization.

## 3.2. Revisiting the Karhunen-Loève decomposition

We assume that the evolution of a certain field that depends on the physical space $\mathbf{x}$ and on time $t$, $u(\mathbf{x},t)$ is known. In practical applications, this field is expressed in a discrete form, that is, it is known at the nodes of a spatial mesh and at some times, i.e. $u(\mathbf{x}_i, t^p) = u_i^p$. We can also write, introducing a spatial interpolation: $u^p(\mathbf{x}) = u(\mathbf{x}, t = p\Delta t)$; $\forall p \in [1, \dots, P]$. The main idea of the Karhunen-Loève (KL) decomposition is how to obtain the most typical or characteristic structure $\varphi(\mathbf{x})$ among these $u^p(\mathbf{x})$, $\forall p$. This is equivalent to obtaining functions $\varphi(\mathbf{x})$ maximizing $\alpha$

$$ \alpha = \frac{\sum_{p=1}^{p=P} \left[ \sum_{i=1}^{i=N} \varphi(\mathbf{x}_i) u^p(\mathbf{x}_i) \right]^2}{\sum_{i=1}^{i=N} (\varphi(\mathbf{x}_i))^2} \quad (14) $$

which leads to:

$$ \sum_{p=1}^{p=P} \left[ \left[ \sum_{i=1}^{i=N} \tilde{\varphi}(\mathbf{x}_i) u^p(\mathbf{x}_i) \right] \left[ \sum_{j=1}^{j=N} \varphi(\mathbf{x}_j) u^p(\mathbf{x}_j) \right] \right] = \alpha \sum_{i=1}^{i=N} \tilde{\varphi}(\mathbf{x}_i) \varphi(\mathbf{x}_i) \quad (15) $$

$\forall \tilde{\varphi}$, which can be rewritten in the form
$$ \sum_{i=1}^{i=N} \left[ \sum_{j=1}^{j=N} \left\{ \sum_{p=1}^{p=P} u^p(\mathbf{x}_i) u^p(\mathbf{x}_j) \varphi(\mathbf{x}_j) \right\} \tilde{\varphi}(\mathbf{x}_i) \right] = \alpha \sum_{i=1}^{i=N} \tilde{\varphi}(\mathbf{x}_i) \varphi(\mathbf{x}_i) \quad (16) $$

Defining vector $\varphi$ such that its $i$-component is $\varphi(\mathbf{x}_i)$, Eq. (16) takes the following matrix form

$$ \tilde{\varphi}^T \mathbf{k} \varphi = \alpha \tilde{\varphi}^T \varphi ; \quad \forall \tilde{\varphi} \Rightarrow \mathbf{k} \varphi = \alpha \varphi \quad (17) $$

where the two points correlation matrix is given by

$$ \mathbf{k}_{ij} = \sum_{p=1}^{p=P} u^p(\mathbf{x}_i) u^p(\mathbf{x}_j) \Leftrightarrow \mathbf{k} = \sum_{p=1}^{p=P} \mathbf{u}^p(\mathbf{u}^p)^T \quad (18) $$

which is symmetric and positive definite. If we define the matrix $Q$ containing the discrete field history:

$$ Q = \begin{pmatrix} u_1^1 & u_1^2 & \cdots & u_1^P \\ u_2^1 & u_2^2 & \cdots & u_2^P \\ \vdots & \vdots & \ddots & \vdots \\ u_N^1 & u_N^2 & \cdots & u_N^P \end{pmatrix} \quad (19) $$

it is easy to verify that the matrix $\mathbf{k}$ in Eq. (18) results in

$$ \mathbf{k} = Q Q^T \quad (20) $$

### 3.3. Reduced modeling

If the evolution of a certain field is known

$$ u(\mathbf{x}_i, t^p) \equiv u_i^p, \quad \forall i \in [1, \dots, N], \quad \forall p \in [1, \dots, P] \quad (21) $$

from some direct simulations, or from experimental measurements, then matrices $Q$ and $\mathbf{k}$ can be computed and the eigenvalue problem given by Eq. (17) can be solved. The solution of Eq. (17) results in $N$ couples eigenvalue-eigenvector. However, in a large number of models involving regular time evolutions of the solution, the magnitude of the eigenvalues decreases very fast, evidencing that the solution evolution can be represented as a linear combination of a reduced number of functions (the eigenvectors related to the highest eigenvalues).
In our numerical applications we consider the eigenvalues ordered as $\alpha_1 > \alpha_2 > \cdots > \alpha_N$. The $n$ eigenvalues belonging to the interval $\alpha_1 > \cdots > \alpha_n$, with $\alpha_n > \alpha_1 \cdot 10^{-8}$ and $\alpha_{n+1} < \alpha_1 \cdot 10^{-8}$ are selected, because their associated eigenvectors are expected to be sufficient to represent accurately the entire solution evolution. In a large variety of models $n \square N$ and moreover $n$ only depends on the regularity of the solution evolution, but neither on the dimension of the physical space (1D, 2D or 3D) nor on the size of the model ($N$).

The reduced approximation basis consist of the $n$ eigenvectors $\varphi_1, \cdots, \varphi_n$, allowing to define the basis transformation matrix $\mathbf{B}$:

$$ \mathbf{B} = (\varphi_1, \varphi_2, \cdots, \varphi_n) \qquad (22) $$

whose size is $N \times n$. Thus, the vector containing the field nodal values $\mathbf{u}$ can be expressed as:

$$ \mathbf{u} = \sum_{i=1}^{n} \varphi_i \cdot \xi_i(t) = \mathbf{B} \cdot \xi(t) \qquad (23) $$

Now, if we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form

$$ \mathbf{A} \mathbf{u}^p = \mathbf{f}^{p-1} \qquad (24) $$

where $\mathbf{f}^{p-1}$ accounts for the solution at the previous time step, taking into account Eq. (23) it reduces to:

$$ \mathbf{A} \mathbf{u}^p = \mathbf{f}^{p-1} \Rightarrow \mathbf{A} \mathbf{B} \xi^p = \mathbf{f}^{p-1} \qquad (25) $$

and multiplying both terms by $\mathbf{B}^T$ it results

$$ \mathbf{B}^T \mathbf{A} \mathbf{B} \xi^p = \mathbf{B}^T \mathbf{f}^{p-1} \qquad (26) $$

which proves that the final system of equations is of low order, i.e. the dimensions of $\mathbf{B}^T \mathbf{A} \mathbf{B}$ are $n \times n$, with $n \square N$, and the dimensions of both $\xi$ and $\mathbf{B}^T \mathbf{f}^{p-1}$ are $n \times 1$.

**Remark 1.** Equation (26) can be also derived introducing the approximation (23) into the Galerkin form of the partial differential equation.

### 3.4. Reduced basis adaptivity
The just described strategy allows for very fast computation of large size models. For example one could solve the full model using some standard discretization technique (finite differences, finite elements, boundary elements, ...) for a small time interval and then define matrix **Q** and **k** allowing compute the reduced approximation basis transformation **B** that leads to the reduced solution procedure illustrated by Eq. (26). Other possibility consists in solving a model in the whole time interval and then extracting the most representative functions that could be used for solving some “similar” models. We come back to this discussion later.

However, in any case, it is not guaranteed that this reduced basis that was built in the first scenario from the solution known within a short time interval, and in the second one for a particular model different to the present one, remains accurate for describing the solution in the entire simulation interval or for any other “similar” model respectively. In the first case it is obvious that during the simulation material properties, boundary conditions ... could change, compromising the validity of the reduced basis. In the second case, the model being different to the one that served to extract the reduced basis nothing guarantees the validity of that reduced approximation basis.

In this manner, if one would compute reduced model solutions and keep the confidence on the related solution, a check of the solution accuracy must be performed and an enrichment strategy must be defined in order to adapt the reduced approximation basis in order to capture the new events present in the solution evolutions which cannot be described accurately from the original reduced approximation basis.

For this purpose, Ryckelynck proposed (Ryckelynck, 2005) to start with a low order approximation basis, using some simple functions (e.g. the initial condition in transient problems) or using the eigenvectors of a “similar” problem previously solved or the ones coming from a full simulation in a short time interval. Now, we compute *S* iterations of the evolution problem using the reduced model (26) without changing the approximation basis. After these *S* iterations, the complete discrete system (25) is constructed, and the residual **R** evaluated:

$$ \mathbf{R} = \mathbf{A}\mathbf{u}^S - \mathbf{f}^{S-1} = \mathbf{A}\mathbf{B}\xi^S - \mathbf{f}^{S-1} \quad (27) $$

If the norm of the residual is small enough, $||\mathbf{R}|| < \varepsilon$, with $\varepsilon$ a threshold value small enough, we can continue for other *S* iterations using the same approximation basis. On the contrary, if the residual norm is too large, $||\mathbf{R}|| \ge \varepsilon$, we need to enrich the approximation basis and compute again the last *S* iterations. This enrichment is built using some Krylov's subspaces, in our case the three first subspaces: $\mathbf{B} \leftarrow (\mathbf{B}, \mathbf{R}, \mathbf{AR}, \mathbf{A}^2\mathbf{R})$.
One could expect the enrichment process increases continuously the size of the reduced approximation basis, but in fact, after reaching the convergence, a Karhunen-Loève decomposition is performed on the whole past time interval in order to extract the significant information as well as to define an orthogonal reduced approximation basis. The interested reader can refer to Ryckelynck et al. (2006) and the references therein for a more detailed and valuable description of the computational algorithm.

### 3.5. Illustrating the applicability of reduced bases

We are considering in this section for the sake of clarity a simple 1D model (the extension to multidimensional models is straightforward) related to the heat transfer equation (we omit the units, all them being expressed in metric system):

$$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \qquad (28) $$

with $\alpha = 0.01$, $t \in ]0,30]$ and $x \in ]0,1[$. The initial condition reads $T(x,t=0)=1$ and the boundary conditions are given by $-\lambda \frac{\partial T}{\partial x}\bigg|_{x=0,t} = q(t)$ and $-\lambda \frac{\partial T}{\partial x}\bigg|_{x=1,t} = 0$ ($\lambda = 0.01$). The boundary source $q(t)$ is prescribed to different values during the simulations that follow.

Equation (28) is discretized by using the implicit finite element method on a mesh that consists of 100 nodes, where a linear approximation is defined in each of the resulting 99 elements. The time step was set to $\Delta t = 0.1$. The resulting discrete system can be written as:

$$ \mathbf{K}\mathbf{T}^p = \mathbf{M}\mathbf{T}^{p-1} + \mathbf{q}^p \qquad (29) $$

where vector $\mathbf{q}^p$ accounts for the boundary heat flux source at each time step $p$.

**Remark 2.** We use the FEM because of the one-dimensionality of the model, but obviously all the results can be extended to any other discretization technique, and in particular to the BEM previously introduced.

First, we consider the solution of the thermal model described above and related to the following boundary heat source:
$$ q(t) = \begin{cases} 1 & 0 < t < 10 \\ 0 & t \ge 10 \end{cases} \qquad (30) $$

The computed solution is depicted in figure 6, where the temperature profiles at times $t_p = p$ ($p = 1, 2, \dots, 30$) are represented. The evolution in the 10 first seconds (heating stage) is depicted in red. In the remaining time interval no more heating sources exist, hence the heat shifts by a conduction mechanism from hotter towards the coldest zones. The profiles within this time interval are represented in blue.
<figure><img src="image_8.png" /><figcaption>Figure 6 - Temperature profiles related to the thermal model with the source term modeled by Eq. (30).</figcaption></figure>

Now, from these 30 profiles we could define matrices **Q** and **k**, which lead to the eigenvalue problem allowing to extract the significant eigenvectors, according to Eq. (4). The resulting eigenvalues are: $\alpha_1 = 1790$, $\alpha_2 = 1.1$, $\alpha_3 = 0.1$ and $\alpha_j < \alpha_1 \times 10^{-8} (4 \le j \le 100)$. This result implies that the whole solution evolution could be accurately represented as a linear combination of the 3 eigenvectors related to the first 3 highest eigenvalues. In order to impose easily the initial condition, it may be possible to add the initial condition to these eigenvectors (even if then the resulting approximation basis is no more orthogonal). Figure 7 depicts the resulting approximation functions, that consist of the 3 eigenfunctions related to the 3 highest eigenvalues and the initial condition, all them normalized, and referred as $\Phi_j$. These functions allow defining matrix **B** and then the reduced model derived from Eq. (29):

$$ \mathbf{B}^T \mathbf{K} \mathbf{B} \xi^p = \mathbf{B}^T (\mathbf{M} \mathbf{B} \xi^{p-1} + \mathbf{q}^p) \quad (31) $$

that only involves 4 degrees of freedom. Even in the case of non linear models and an implicit discretization, only the inversion of a matrix of size 4 is required at each time step.
If we assume that the initial condition has been placed in the first column of **B**, then the initial condition in the reduced basis writes: $(\xi^0)^T = (1 \ 0 \ 0 \ 0)$. Now, from this condition, Eq. (31) can be applied to compute the whole time evolution. Obviously, the global solution can be obtained from the reduced one according to the basis transformation relationship: $T^p = \mathbf{B}\xi^p$. Figure 8 compares a few temperature profiles obtained using the global model (Eq. (29)) and that were depicted in figure 6, with those obtained using the reduced model (Eq.(31)). An excellent accuracy can be noticed. This accuracy is not surprising because, as indicated before, the four approximation functions used were the ones related to the highest eigenvalues and that, consequently represent the optimal reduced approximation basis.

<figure><img src="image_9.png" /><figcaption>**Figure 7 -** Functions defining the reduced approximation basis.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8** - Complete (continuous line) versus reduced (stars) model solutions.</figcaption></figure>

In order to conclude about the applicability of this reduced approximation basis to simulate problems different from the one that served to compute it, we consider the thermal model defined in the same domain, with the same initial condition, but with a slightly different boundary heat source term:

$$q(t) = \begin{cases} \frac{t}{20} & 0 < t < 20 \\ \frac{t-30}{5} & t \ge 20 \end{cases} \qquad (32)$$

Figure 9 compares the reference solution (continuous line) computed using the model represented by Eq. (29) with that obtained using the reduced model (31) – stars- but where the reduced approximation basis consists of the four functions represented in Fig. 7 and that were associated with the thermal model related to the boundary condition given by Eq. (30). We can notice the excellent accuracy, somewhat unexpected, because the non evident compatibility between the problem solutions defined by Eqs. (30) and (32), and then the unexpected ability of the approximation functions extracted from the solution of the thermal model defined by Eq. (30) to describe the solution of the thermal model related to Eq. (32).
<figure><img src="image_11.png" /><figcaption>**Figure 9** - Complete (continuous line) versus reduced (stars) model solutions related to the thermal model associated to the thermal source given by Eq. (19).</figcaption></figure>

From this result we can start to realize the potentiality of model reduction, but in any case two questions remain open: (i) how to quantify the quality of a reduced solution without the necessity of computing the global solution?; and (ii) in the case of noticing a lack of accuracy, how to enrich the reduced approximation basis in order to improve the solution accuracy?

To address these questions, we consider the technique originally proposed by Ryckelynck (Ryckelynck, 2005) that consists in computing the solution residual defined at a certain time step by

$$R = KT^p - MT^{p-1} - q^p \qquad (33)$$

This residual can be used to quantify the accuracy of the reduced solution, and allows addressing the first question. Now, concerning the second one, we are assuming that the residual resulting from the application of Eq. (20) is greater than a threshold value, i.e. $||R|| \ge \varepsilon$. A natural choice consists in enriching the reduced approximation basis by adding this residual (that is orthogonal in a Galerkin sense to the approximation functions) and some of the Krylov's subspaces related to it (according to the procedure previously described).
For the thermal model related to the boundary source given by Eq. (32) and the reduced approximation basis depicted in Fig. 7, the residual norm at $t = 30$ was $\|\mathbf{R}\| = 2.6 \times 10^{-5}$, justifying the excellent agreement between the global and reduced solutions noticed in Fig. 9.

Figure 10 depicts the normalized residual $(\mathbf{R} \leftarrow \mathbf{R}/\|\mathbf{R}\|)$. Even if the residual is very small, it can be noticed that the largest deviations are concentrated around the left boundary, where the boundary thermal source applies.

<figure><img src="image_12.png" /><figcaption>**Figure 10** - Normalized residual at $t = 30$ related to the solution associated with Eq. (32) and computed from the reduced approximation basis depicted in Fig. 7.</figcaption></figure>

Now, we proceed to enrich the reduced approximation basis by introducing this residual (the first Krylov's subspace) into matrix **B** according to: $\mathbf{B} \leftarrow (\mathbf{B}, \mathbf{R})$. If the thermal model related to Eq. (32) is solved again, but using the just updated (enriched) reduced approximation basis, then the norm of the residual decreases, justifying the introduction of the Krylov's subspaces related to the residual to improve the reduced solution accuracy.

## 3.6. Discussion

Similar results could be obtained by considering more complex thermal models. We have analyzed several different scenarios. One of them consisted of a thermal model
that was solved for a particular evolution of the volumetric source term (the boundary conditions were assumed fixed in this case). A first solution allowed us to define the reduced approximation basis by solving the eigenvalue problem and selecting the eigenvectors related to the highest eigenvalues. Now, different evolutions in time and space of the volumetric source term were considered and the associated solutions computed by using the same reduced basis. The computed solutions were compared with the reference ones obtained by using the complete approximation basis. In all the analyzed cases the agreement was very good, always improved by adding to the reduced basis the computed residuals associated to the reduced basis solutions, as just described.

The previous numerical examples illustrate the model reduction capabilities. The main originality lies in the capability of checking for the solution accuracy and the eventual possibility to adapt the reduced basis by adding the residual and some Krylov's subspaces generated by it into the reduced approximation basis. These appealing capabilities were exploited to solve many models (see the works referenced in the introduction to this section). We have also experienced that models involving weak fixed discontinuities (thermal models in non-homogeneous media) accept reduced approximation basis, from which the evolving solutions can be accurately represented. When these discontinuities are evolving within the domain, the situation becomes a bit more complicated. If general, strong discontinuities moving within the domain don't accept a reduced description. That is, if we know the solution evolution and we apply the Karhunen-Loève decomposition, most of the eigenvalues must be retained in the model description. This situation seems quite obvious. Other difficulties are found when the models imply moving meshes (updated Lagrangian formulations). In that case, a reduced basis can be extracted (Ryckelynk et al., 2006), but it cannot be used to solve models in which the evolution of the nodal positions differ from the one that served to define the reduced basis. The fact that level-set based descriptions of moving interfaces accept reduced approximations opens new perspectives within the framework of the partition of unity based discretizations.

In what follows we focus on a direct consequence of the examples discussed in the previous section, concerning the optimization. As it is well known, in order to perform optimization one needs a minimization strategy (first or higher order) and a direct solver that is called for each tentative state of the design parameters. Because one must solve numerous direct problems, the computing time of such direct solver becomes crucial. The simplest alternative lies in solving the complete model for one (or some) point within the design space and extract the reduced basis, and then for any other point (given by the minimization strategy) computed the solution from that reduced approximation basis. In any case, the solution can be improved by enriching the reduced bases by applying the strategy described above. In the next sections we illustrate the capabilities of such one procedure.
# 4. Application of the DRBEM reduced model to mold cooling optimization

## 4.1. Reduced model coupled with DRBEM

We solve the eigenvalue problem defined in section 3 selecting the eigenfunctions $\phi_n$ associated with the eigenvalues $\alpha_n (n \in [1, N])$ belonging to the interval defined by the largest eigenvalue, such as $\phi_n$'s sum is upper or equal to 99.9% of $\phi_N$'s sum (N number of nodes). In practice, n is much lower than N. The **B** matrix is then assembled and used to approximate the temperature. The main steps of the direct simulation coupled with the reduced model are summarized on figure 11.

<figure><img src="image_13.png" /><figcaption>Figure 11 – Direct DRBEM simulation with reduced model</figcaption></figure>

## 4.2. Overall optimization methodology

We present in this section how we formulate the problem under a mathematical programming form. In the sequel, $x$ will denote the vector of optimization variables (position and shape parameters for the cooling channels). Since the output of the heat-transfer problem is a function of $x$, we shall make the dependence explicit of the temperature measurements upon the position and shape parameters. Most practical optimization problems involve several (often contradictory) objective functions. The simplest way to proceed in such a multi-criterion context is to consider as objective function a weighted sum of the various criteria. This involves choosing appropriate weighting parameter values. An obvious alternative is to use one criterion as objective function while requiring, in the constraints, maximal threshold levels for the remaining criteria. We choose here the latter approach because we do know a
threshold level value for the maximal temperature variation under which any variation is equally acceptable. More precisely, we formulate our problem under the form:

$$ \min_{x} \|T(x)\| \quad \text{subject to } f(T(x)) \text{ and } g(x) \le 0 \qquad (34) $$

where f is a real-valued function used to stipulate the uniformity-temperature constraint, and g(x) is a general vector-valued non-linear function. The complete methodology to couple the thermal solver and the optimization algorithm procedure is presented in (Pirc, 2009).

The general constraints g(x) ≤ 0 represent any geometry related or other industrial constraints, such as:

* upper/lower-bound constraints on the xᵢ,

* keeping the cooling channels within the mold,

* technically-forbidden zones where we cannot position the cooling channels (for instance due to the presence of ejectors),

* constraints stipulating a minimal distance between every pair of cooling channels to avoid inter-channels collision.

### 4.3. Application

In this section, we report numerical simulations on a 3D plastic part whose features are displayed on Figure 12 (units in mm). It is a semi-industrial injection mold design for the European project: Eurotooling 21.

<figure><img src="image_14.png" /><figcaption>Figure 12 – Plastic part dimensions</figcaption></figure>

The mold is meshed using 5592 linear triangles and each cooling channel using 340 quadrangles (figure 13).
<figure><img src="image_15.png" /><figcaption>Figure 13 – Upper-part of the mold mesh</figcaption></figure>

The thermo-physical properties of the polymer as well the mold material are referenced in table 1. The boundary conditions are the same as defined in section 2.2.

<table><caption>Table 1 – Thermo-physical properties</caption><thead><tr><th></th><th>Polymer (PP)</th><th>Mold (Steel)</th></tr></thead><tbody><tr><td>λ [W.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>0.63</td><td>34</td></tr><tr><td>ρ [kg.m<sup>-1</sup>]</td><td>891</td><td>7800</td></tr><tr><td>Cp [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>2740</td><td>460</td></tr></tbody></table>

The history matrix, corresponding to the first injection cycle time, is computed using transient DRBEM code. The mold temperature, for the next injection cycle time, is computed using the reduced model. The optimization objective consists in minimizing the maximal temperature while minimizing temperature variations:

$$ \text{minimize } \max_{i \in N}(T_i) \quad \text{subject to} \quad \sum_{i \in N} |T_i - T_{av}| \le \sigma \qquad (35) $$

where N is the number of elements and $T_{av}$ the average surface temperature. For illustration purposes, we consider here 8 cooling channels and the constraints and optimization variables are sketched in figure 14.
<figure><img src="image_16.png" /><figcaption>**Figure 14** –Constraints and optimization variables</figcaption></figure>

The geometrical optimization parameters are the coordinates of the end points, $P_1$ and $P_2$ of each cooling channel (figure 14). Since $P_2$ can be expressed in terms of the other coordinates and since the channel length ($L$) is constant, the optimization parameters for locating the $i$-th cooling channel are completely determined by $P_i^i = (X_i, Y_i, Z_i) \quad i=1,...,8$. For this application, $Z_i$ is fixed and therefore the problem reduces to 16 optimization variables.

We use as starting point a heuristic solution provided by an experienced engineer. On average, one objective function evaluation requires 14 min of CPU time (one direct computation). Since we compute gradients using finite difference approximation (associated with SQP method), one optimization iteration involves 4 hours of CPU time on *Macintosh 1.83 GHz Intel Core 2 Duo*. 24 optimization iterations were necessary in order to achieve convergence for one injection cycle (figure 15).

<figure><img src="image_17.png" /><figcaption>**Figure 15** –Convergence history</figcaption></figure>
In addition, the surface temperature distribution of the mold is presented in figure 16, with the temperature profile at the mold surface before and after optimization are shown in figure 17. We observe that both temperature variance and temperature average decrease significantly.

<figure><img src="image_18.png" /><figcaption>Figure 16 –Surface temperature distribution of the mold</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17 –Temperature profile at the surface of mold cavity before and after optimization (z=0.02 m)</figcaption></figure>

We need 96 hours to perform a complete optimization without the reduction model. If we use now the reduction model and DRBEM, we reduce CPU time to 7h40 (one direct computation of 14 minutes and 24 optimization iterations of 18 minutes approximately). CPU time is divided approximately by 13.

**4.4. Conclusion**

We introduced a methodology based on the use of DRBEM to solve the 3D heat transfer equation during the cooling step of the injection molding process. The preliminary computation tests on a semi-industrial plastic part showed that the approach is viable for optimizing the design of cooling channels for injection molding. The numerical modeling and optimization methodology can easily take into account a large range of industrial constraints. Various optimization criteria can be provided by the user (either directly as a cost function or within constraints).
Another interesting aspect consists in using this technique in order to compute multi-cycles injection mold cooling. The reduction model allows to extract from the first cycle the relevant eigenfunctions associated with the eigenvalues and consequently to calculate very rapidly all the other cycles.

## REFERENCES

J.-F. Agassant, P. Avenas, J.-Ph. Sergent, P.J. Carreau. *Polymer Processing – Principle and Modelling*, Hanser Publishers, New York 1991

A. Ammar, D. Ryckelynck, F. Chinesta, R. Keunings. *On the reduction of kinetic theory models related to finitely extensible dumbbells*. *Journal of Non-Newtonian Fluid Mechanics*, **134**, 136-147, 2006.

A. Ammar, E. Pruliere, J. Ferec, F. Chinesta, E. Cueto. Coupling finite elements and reduced approximation bases. *European Journal of Computational Mechanics*. In press, 2009.

A. Bialecki, P. Jurgas, G. Kuhn. *Dual reciprocity BEM without matrix inversion for transient heat conduction*. *Engineering Analysis with Boundary Element*, **26**, 227-236, 2002.

A. Bikas, A. Kanarachos. The dependence of cooling channels system geometry parameters on product quality as a result of uniform mold cooling. *Proceedings of the 57th Annual Technical Conference ANTEC 99*, 578–583, 1999.

E. Boillat, R. Glardon, D. Paraschivescu. *Journal de Physique IV*, **102**, 27-38, 2002.

C.A. Brebbia, J. Dominguez. Boundary elements, an introductory course. WIT Press/Computational Mechanics Publications, 1992.

C. A. Brebbia, C. S. Chen, H. Power. *Dual reciprocity method using compactly supported radial basis functions*. *Communications in Numerical Methods in Engineering*, **15/2**, 137-150, 1999.

F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene, F. Boust. Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization. *Computer Methods in Applied Mechanics and Engineering*, **197/5**, 400-413, 2008.

R. Fletcher. *Practical methods of optimization*, Wiley, 1987.

L. Godinho, A. Tadeu, N. Simoes. Study of transient heat conduction in 2.5D domains using the boundary element method. *Engineering Analysis with Boundary Elements*, **28/6**, 593–606, 2004.
P.J. Holmes, J.L. Lumley, G. Berkoozld, J.C. Mattingly, R.W. Wittenberg. Low-dimensional models of coherent structures in turbulence. Physics Reports, **287**, 1997.

J. Huang, G. Fadel, Journal of Mechanical Design, **123**, 226–239, 2001.

K. Karhunen. Über lineare Methoden in der wahrscheinlichkeitsrechnung. Ann. Acad. Sci. Fennicae, ser. Al. Math. Phys., **37**, 1946.

P. Krysl, S. Lall, J.E. Marsden. Dimensional model reduction in non-linear finite element dynamics of solids and structures. Int. J. Numer. Meth. in Engng., **51**, 479–504, 2001.

M.M. Loève. Probability theory. The University Series in Higher Mathematics, 3rd Ed. Van Nostrand, Princeton, NJ, 1963.

E.N. Lorenz. Empirical orthogonal functions and statistical weather prediction. MIT, Department of Meteorology, Scientific Report N1, Statistical Forecasting Project, 1956.

Y. Lu, D. Li, J. Xiao. Simulation of cooling process of injection molding. Progress in Natural Science, **6/2**, 227–234, 1996.

E. Mathey, L. Penazzi, F.M. Schmidt, F. Ronde-Oustau. Automatic optimization of the cooling of injection mold based on the boundary element method. Materials Processing and Design: Modeling, Simulation and Applications. Proceedings NUMIFORM 04, 222–227, 2004.

J. Moller, M. Carlson, R. Alterovitz, J. Swartz. Post-ejection cooling behavior of injection molded parts. Proceeding of the 56th Annual Technical Conference ANTEC 98, 525–527, 1998.

S. Niromandi, I. Alfaro, E. Cueto, F. Chinesta. Real-time deformable models of non-linear tissues by model reduction techniques. Computer Methods and Programs in Biomedicine, **91**, 223–231, 2008.

S.W. Opolski, T.W. Kwon. Injection molding cooling system design. Annual Technical Conference and Exhibition, Society of Plastics Engineers, 264–268, 1987.

T.A. Osswald. Polymer Processing – Fundamentals. Hanser Publishers, 1998.
H.M. Park, D.H. Cho. The Use of the Karhunen-Loève decomposition for the modelling of distributed parameter systems, Chem. Engineer. Science, **51**, 81-98 (1996).

S. Park, T. Kwon, Polymer Engineering and Science, **38**, 1450–1462, 1998.

R. Pasquetti, D. Petit. Inverse diffusion by boundary elements. Engineering Analysis with Boundary Elements, **15**, 197–205, 1995.

N. Pirc, F. Bugarin, F.M. Schmidt, M. Mongeau. 3D BEM-based cooling-channel shape optimization for injection molding processes. International Journal for Simulation and Multidisciplinary Design Optimization, 2/3, 245-252, 2008.

N. Pirc, F.M. Schmidt, M. Mongeau, F. Bugarin, F. Chinesta. Optimization of 3D Cooling Channels in injection molding using DRBEM and Model reduction. 12th ESAFORM Conference on Material Forming, 2009.

N.S. Rao, G. Schumacher, N.R. Schott, K.T. O’Brien. Optimization of cooling systems in injection molds by an easily applicable analytical model. Journal of Reinforced Plastics and Composites, **21/5**, 451–459, 2002.

D. Ryckelynck. A priori hyperreduction method: an adaptive approach. Journal of Computational Physics, **202**, 346-366, 2005.

D. Ryckelynck, L. Hermanns, F. Chinesta, E. Alarcón. An efficient “a priori” model reduction for boundary element models. Engineering Analysis with Boundary Elements, **29**, 796-801, 2005.

D. Ryckelynck, F. Chinesta, E. Cueto, A. Ammar. On the “a priori” model reduction: Overview and recent developments. Archives of Computational Methods in Engineering, State of the Art Reviews, **13/1**, 91-128, 2006.

L. Sirovich. Turbulence and the dynamics of coherent structures part I: Coherent structures. Quaterly of applied mathematics, XLV, 561–57, 1987.

Sutradhar, G.H. Paulino, L.J. Gray. Transient heat conduction in homogeneous and nonhomogeneous materials by the Laplace transform Galerkin boundary element method. Engineering Analysis with Boundary Elements, **26/2**, 119–132, 2002.

L. Tang, K. Pochiraju, C. Chassapis, S. Manoochehri. Journal of Mechanical Design, **120**, 165–174, 1998.
N. Verdon, C. Allery, C. Béghein, A. Hamdouni, D. Ryckelynck. Reduced-Order Modelling for solving linear and non-linear equations. Communications in Numerical Methods in Engineering. In press, 2009.

M. H. Wesselmann, Impact of moulding conditions on the properties of short fibre reinforced high performance thermoplastic parts. PhD Thesis, Ecole des Mines Albi, 1998.

S.Y. Yang, L. Lien. Effects of cooling time and mold temperature on quality of moldings with precision contour. Advances in Polymer Technology, 15/4, 289-295, 1996.