# A Simulation App based on reduced order modeling for manufacturing optimization of composite outlet guide vanes

Jose Aguado, Domenico Borzacchiello, Chady Ghnatios, François Lebel, Ram Upadhyay, Christophe Binetruy, Francisco Chinesta

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

## To cite this version:

Jose Aguado, Domenico Borzacchiello, Chady Ghnatios, François Lebel, Ram Upadhyay, et al. A Simulation App based on reduced order modeling for manufacturing optimization of composite outlet guide vanes. Advanced Modeling and Simulation in Engineering Sciences, 2017, 10.1186/s40323-017-0087-y. [https://hal.science/hal-01925863](https://hal.science/hal-01925863)

HAL Id: hal-01925863

https://hal.science/hal-01925863

Submitted on 18 Nov 2018

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
# A *Simulation App* based on reduced order modeling for manufacturing optimization of composite outlet guide vanes

José V. Aguado

Domenico Borzacchiello

Chady Ghnatios

François Lebel

Ram Upadhyay

Christophe Binetruy

Francisco Chinesta

## Abstract

Composites manufacturing processes usually involve multiscale models in both space and time, highly non-linear and anisotropic behaviors, strongly coupled multiphysics and complex geometries. In this framework, the use of simulation for real-time decision making directly in the manufacturing facility is still precluded nowadays, in spite of the impressive progresses reached in numerical analysis and computer science during the last decade. In this paper, a process-specific simulation tool based on reduced order modeling is introduced, the *Simulation App*. This concept is presented through a practical case involving a multi-physics and coupled problem describing the manufacturing process of a composite outlet guide vane. We show that several manufacturing settings can be simulated in few seconds with the *Simulation App*, thus enabling fast process optimization. Finally, the advantages over general-purpose simulation software, in the context of process simulation, are discussed.

Simulation App; Reduced order modeling; Real time simulation; Numerical simulation; Composites forming processes; Curing; Consolidation; Proper Generalized Decomposition; Parametric solutions; Computational Vademecums

## 1 Introduction

Efficient simulation of composite manufacturing processes remains even nowadays, in many cases, a challenging issue, mainly when they involve rich 3D behavior, multi-physics and the necessity of solving many scenarios very fast for optimization purposes [1, 2]. Composites manufacturing processes involve many different physics. For instance, impregnation of fibrous reinforcements involves flow models through porous media, combined with the mould compression for ensuring an appropriate degree of consolidation [3, 4]. During consolidation the resin curing kinetics strongly affects its viscosity and therefore the flow conditions [5].
The intimate coupling between the resin flow, the reinforcement permeabil-
ity, that decreases upon compression of the mould, and the viscosity, that de-
pends on the temperature and the curing degree, can be at the origin of dif-
ferent process defects. Large viscosities imply high pressures that can generate
irreversible preform deformations and even local reinforcement displacements
creating local inhomogeneity defects. Moreover, residual stresses are another
manifestation of this intimate multi-physics coupling [1, 6].

Another important issue encountered in the simulation of composites man-
ufacturing is the one related to the process control and optimization [7, 8, 9].
In general, optimization implies the definition of a cost function and the search
of the optimum process parameters defining the minimum of that cost function.
The procedure starts by a guessed set of process parameters. Then the process
is simulated using suitable numerical methods. The solution of the model is the
most expensive step of the optimization procedure. As soon as the solution is
available, the cost function can be evaluated and its optimality checked. If the
chosen parameters do not define a minimum (at least local) of the cost func-
tion, the process parameters should be updated and the solution recomputed.
The procedure continues until reaching the minimum of the cost function. The
solution of the process model is a tricky task that demands important compu-
tational resources and usually implies extremely large computing times. Thus,
usual optimization procedures are inapplicable under the real-time constraint.
The same issues are encountered when dealing with inverse analysis in which
material or process parameters are expected to be identified from numerical sim-
ulation, by looking for the unknown parameters so that computed fields match
the ones measured experimentally.

Until now, the solution consisted in using the more and more powerful com-
puting platforms and techniques for speeding up standard discretization tech-
niques. Appealing alternatives for circumventing, or at least alleviating, these
issues lie in the use of reduced order modeling (ROM) [10, 11, 12] strategies.
ROM is based on the observation that the family of parametric solutions of a
given model usually contains much less information than it was originally as-
sumed when the discrete model was built. Proper Orthogonal Decomposition,
Proper Generalized Decomposition and Reduced Basis are nowadays widely
considered from a fundamental and applicative viewpoints.

Proper Orthogonal Decomposition (POD) is a general technique for extract-
ing the most significant characteristics of a system’s behavior and representing
them in a set of “POD basis vectors” [13, 14]. These basis vectors then provide
an efficient (typically low-dimensional) representation of the key system behav-
ior, which proves useful in a variety of ways. The most common use is to project
the system governing equations onto the reduced-order subspace defined by the
POD basis vectors. This yields an explicit POD reduced model that can be
solved in place of the original system. The POD basis can also provide a low-
dimensional description in which to perform parametric interpolation [15, 16],
infill missing or “gappy” data [19], hyper-reduced approximations [17, 18], and
perform model adaptation. There is an extensive literature on POD showing it
has broad application across fields. Some review of POD and its applications
can be found in [20, 21] and the references therein.

Another family of ROM techniques is based on the use of a reduced basis
constructed by combining a greedy algorithm and “a posteriori” error indicators
[22, 23, 24]. As for the POD, the Reduced Basis method requires some amount
offline work, but then the reduced basis model can be used online for solving
different models with control of the solution accuracy, because the availability
of error bounds[25]. When the error is unacceptably high, the reduced basis can
be enriched by invoking a greedy adaption strategy. Useful review works on the
subject are [26, 27].

Techniques based on the use of separated representations are at the heart of
the so-called Proper Generalized Decomposition (PGD) methods [28, 29, 30].
Such separated representations were considered in computational mechanics for
separating space and time in transient solutions [31, 32]. Separated representa-
tions were employed for solving multidimensional models suffering the so-called
curse of dimensionality [33, 34, 35] and in the context of stochastic modeling
[36]. Then, they were extended for separating space coordinates making possi-
ble the solution of models defined in degenerated domains, e.g. plate and shells
[37] as well as for addressing parametric models where model parameters were
considered as model extra-coordinates, making possible the offline calculation of
the parametric solution that can be viewed as a metamodel or a computational
vademecum, to be used online for real time simulation, optimization, inverse
analysis and simulation-based control [38]. Some applications in the context
of composites manufacturing processes were addressed in [39] while the multi
physics coupling was successfully achieved in [40].

This work is intended to propose a first approach using reduced order mod-
eling techniques to enhance adaptability of composite manufacturing process
to changeable material and process environments through increased parametric
modeling capabilities. The process selected is the consolidation and curing of a
real part. Consolidation and curing of thermoset pre-impregnated fibers (from
now on, simply referred as pre-pregs) involve different physics: heat transfer,
compression of fiber beds, resin flow and chemical reaction. Strong couplings
exist between these physics and many material parameters come into play. In
this paper we combine several different modeling and simulation strategies for
the efficient solution of a generic multi-physic and coupled problem. In partic-
ular, we propose a ROM-based segregated approach (rather than a monolithic
one) for treating each physics separately. This approach allows applying, in each
case, the most convenient ROM technique. Then, coupling is made by defining
an appropriate parametrization of the coupling variables. A strategy for cou-
pling ROMs of different kind (i.e. reduced solvers and parametric solutions)
is also proposed. The integration of all these models constitutes a *Simulation
App* that allows real-time evaluation of any process conditions. The described
methodology can be extended and generalized to other processing technologies.

The rest of the paper is organized as follows. Section 2 presents the man-
ufacturing process as well as the equations describing the physics involved in
the manufacturing process of composite outlet guide vanes (OGV). Section 3 de-
scribes the simulation strategy, including the different simulation modules based
<figure><img src="image_3.png" /><figcaption>Figure 1: Global view of the OGV part geometry with cross sectional slice (left) and detail of the computational mesh (right).</figcaption></figure>

on ROM and their coupling. In particular, we recall the main features of the ROM methods that were implemented, emphasizing both the construction of the reduced models in the offline stage and their utilization in the online stage, as a part of the *Simulation App*. Finally, Section 4 presents how the simulation modules were integrated into a process-specific simulation tool, a *Simulation App*. We describe the different functionalities of the application, organized by tabs, allowing to the user to perform the usual pre-processing, simulation and visualization operations. The advantages of process-specific simulation tools, such as the *Simulation App*, over general-purpose simulation software are discussed, especially in process optimization framework.

# 2 Process description and physics modeling

The concept of a *Simulation App* is presented in this paper through a practical application: the manufacturing process of a composite blade, and more specifically, an outlet guide vane. In this section we are first describing the manufacturing process from a technological point of view. Then we shall present the physics modeling and constitutive behavior for each considered physical phenomena.

## 2.1 Process description

Figure 1 shows the geometry and the computational mesh of the composite OGV part, which is manufactured by press forming and curing of 66 thermoset unidirectional continuous fiber pre-preg layers, i.e. fibers already impregnated with resin in a partially cured stage.

During the forming process, the composite lay-up is heated by conduction from the top and bottom metallic mold walls and consolidated under press. See Figure 2 for details. The heat initiates the cure reaction and the applied pressure provides the force needed to drain the excess of resin out of the composite, to consolidate pre-preg layers and to reduce voids by compressing the air inside. The temperature raise determines the onset of the cure reaction. Because of the exothermic effects of the curing process and the variability of the kinetic
<figure><img src="image_4.png" /><figcaption>Figure 2: Global view of the press (left) and detail of an OGV part being produced (right). Images courtesy of General Electric Global Research.</figcaption></figure>

and thermal properties of the resin with the temperature and curing degree, the process is highly nonlinear. Moreover, the resin undergoes a strong rheological modification because its viscosity also depends on the temperature and the degree of cure. Therefore the flow conditions vary continuously with time.

The nominal heating cycle is designed so as to apply a constant temperature of 438 K (330 F) on both top and bottom mold walls. The press closure is initially based on a constant closure rate of 0.1524 mm/min (6 mils/min). As the cure reaction advances, the viscosity increases and so the closing force to be applied by the press increases. When the maximum closing force of the press is attained, the control switches to a force-based one. The closure rate is therefore adapted so as to keep the closing force constant. A final technological restriction to be taken into account is the so-called *hard-stop condition*, occurring when the cumulated displacement reaches its maximum value, 3.3 mm (130 mils).

**Remark 1** (On the simulation of the press control system). *Observe that the manufacturing process simulation will be nonlinear due to the just described press control system. If a force-based control applies, one has to solve a nonlinear problem in order to find out the closure rate that keeps the force constant at its maximum value.*

## 2.2 Thermo-kinetic model

The model governing thermo-kinetic is given by the following partial differential equations:

$$ \rho C_p(T, \alpha) \frac{\partial T}{\partial t} = \nabla \cdot \lambda \nabla T + \rho \Delta H \dot{\alpha} \quad (1a) $$

$$ \dot{\alpha} = f^{\text{kin}}(T, \alpha), \quad (1b) $$

where $T(\mathbf{x}, t)$ and $\alpha(\mathbf{x}, t)$ represent the temperature and curing degree fields, respectively. We shall refer to Equation (1a) and to Equation (1b) as the heat and kinetic equations, respectively. The curing function, denoted by $f^{\text{kin}}$, is
given by the phenomenological model by Kamal and Sourour [41], widely used to describe the conversion kinetics of epoxy systems:

$$f^{\text{kin}}(T, \alpha) = \left( k_1 \exp\left(\frac{E_1}{RT}\right) + k_2 \alpha^m \exp\left(\frac{E_2}{RT}\right) \right) (1 - \alpha)^n, \quad (2)$$

Heat convection is neglected because of the creeping flow approximation and it is assumed there is no dispersion effect, i.e. resin and fiber share the same temperature. The coefficients $\rho$ and $C_p$ appearing in Equation (1a) denote the density and the specific heat, while $\lambda$ stands for the thermal conductivity tensor of the composite part. These were obtained from the properties of the resin and fibers using a linear mixing rule. For the specific heat this reads:

$$C_p(T, \alpha) = V_f C_{pf} + (1 - V_f) C_{pr}, \quad (3)$$

where $C_{pf}$ and $C_{pr}$ are the specific heat of the fibers (depending on the temperature) and the resin, respectively. At the initial conditions, i.e. $V_f = V_{f0} = 57\%$, it can be rewritten as follows:

$$C_p(T, \alpha) = (p_1 T + p_0) \alpha + (q_1 T + q_0) (1 - \alpha), \quad (4)$$

where coefficients $p_0$, $p_1$, $q_0$ and $q_1$ are defined in Table 1. Note that only this quantity is varying with the temperature, whereas $\rho$ and $\lambda$ can be reasonably approximated as constants in the system considered in this application. The thermal properties, as well as the coefficients of the kinetic model, were obtained from experimental characterization and are summarized in Table 1.

The solution of the thermo-kinetic model requires defining appropriate initial and boundary conditions. The initial conditions are given by $T(\mathbf{x}, t = 0) = T_0$ and $\alpha(\mathbf{x}, t = 0) = \alpha_0$. It is important to note that the heat equation, Equation (1a), is global in space whereas the kinetic equation, Equation (1b) is local, that is, the solution at a particular location of the part only depends on the thermo-kinetic history at that location. Therefore, boundary conditions are only needed for the heat conduction equation, and in particular they concern the temperature cycle prescribed at the top and bottom surfaces $S^+$ and $S^-$, see Figure 3, that read $T(\mathbf{x} \in S^+, t) = T_t(t)$ and $T(\mathbf{x} \in S^-, t) = T_b(t)$, respectively. The heat losses through the lateral surfaces $\mathcal{L}$ (see Figure 3) can be neglected, that results in $\nabla T(\mathbf{x}, t) \cdot \mathbf{n} = 0$, being $\mathbf{n}$ the unit outwards vector defined on $\mathcal{L}$.

## 2.3 Consolidation model

During the consolidation under press, the excess of resin is drained out of the composite. Assuming that only the resin moves between fibers, that is, preform is compressed but remains nearly at rest, the resin flow model can be described from the Darcy's flow model:

$$\left\{ \begin{array}{l} \nabla \cdot \mathbf{v} = 0 \\ \mathbf{v} = -\eta^{-1} \mathbf{K} \cdot \nabla P \end{array} \right. , \quad (5)$$

where $\mathbf{v}$, $P$ and $\eta$ are the velocity, pressure and viscosity fields and $\mathbf{K}$ is the permeability tensor. The solution of the flow model requires adequate boundary
<figure><img src="image_5.png" /><figcaption>Figure 3: Boundary conditions description for the thermo-kinetic problem.</figcaption></figure>

conditions. In particular it is assumed that in the lateral boundaries $\mathcal{L}$, see Figure 3, the pressure vanishes. On the other hand, it is assumed that the velocity at the top surface corresponds to the compression rate $\dot{U}$ imposed by the press, which is assumed to be directed vertically, while the resin velocity vanishes at the bottom surface. In summary:

$$ \left\{ \begin{array}{l} P(\mathbf{x} \in \mathcal{L}) = 0 \\ \mathbf{v}(\mathbf{x} \in \mathcal{S}^{+}) = (0, 0, -\dot{U})^{T} . \\ \mathbf{v}(\mathbf{x} \in \mathcal{S}^{-}) = \mathbf{0} \end{array} \right. \qquad (6) $$

As the curing reaction advances, the resin viscosity $\eta$ tends to increase until the flow is no longer possible or will induce fiber washing. This behavior is taken into account in the following chemo-rheological model:

$$ \eta(T, \alpha) = \eta_0 \exp\left(\frac{B}{T}\right) \exp(C\alpha^d). \qquad (7) $$

The viscosity evolution is therefore the main link between Equation (1) and Equation (5): the thermo-kinetic simulation influences the pressure field through the viscosity, which depends on both the temperature and the curing degree. Note also that, strictly speaking, the thermo-kinetic simulation depends on the consolidation simulation, because, as the excess of resin is drained out under the action of the press, the fiber volume fraction changes, and so the specific heat does, see Equation (3). However, this link will be neglected in practice thanks to the semi-implicit linearization scheme that was implemented, see Section 3.4 for details.

Finally, the permeability is also assumed to depend on the fiber volume fraction evolution. The principal components obey the following expressions[42]:

$$ K_{11} = K_{22} = c_{\parallel} \frac{(1 - V_f)^3}{V_f^2} \quad \text{and} \quad K_{33} = c_{\perp} \left( \sqrt{\frac{V_{fa}}{V_f}} - 1 \right)^{5/2} . \qquad (8) $$
<table><tr><td>k<sub>1</sub></td><td>1.77 &middot; 10<sup>-5</sup> s<sup>-1</sup></td></tr><tr><td>k<sub>2</sub></td><td>9.30 &middot; 10<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>E<sub>1</sub>/R</td><td>9,200 K</td></tr><tr><td>E<sub>2</sub>/R</td><td>6,450 K</td></tr><tr><td>n</td><td>0.91</td></tr><tr><td>m</td><td>0.37</td></tr><tr><td>&Delta;H</td><td>570 J/g</td></tr><tr><td>&eta;<sub>0</sub></td><td>3.10 &middot; 10<sup>-12</sup> Pa s</td></tr><tr><td>B</td><td>10800 K</td></tr><tr><td>C</td><td>27</td></tr><tr><td>d</td><td>1.77</td></tr><tr><td>c<sub>&parallel;</sub></td><td>9.43 &middot; 10<sup>-13</sup> m<sup>2</sup></td></tr><tr><td>c<sub>&perp;</sub></td><td>2.50 &middot; 10<sup>-12</sup> m<sup>2</sup></td></tr><tr><td>V<sub>fa</sub></td><td>78.5%</td></tr></table>

<table><caption>Table 1: Material parameters considered in the physical modeling of the OGV part manufacturing process. Courtesy of GE Global Research.</caption><tr><td>p<sub>0</sub></td><td>0.1203 J/gK</td></tr><tr><td>p<sub>1</sub></td><td>0.005 J/gK<sup>2</sup></td></tr><tr><td>q<sub>0</sub></td><td>0.5905 J/gK</td></tr><tr><td>q<sub>1</sub></td><td>0.0019 J/gK<sup>2</sup></td></tr><tr><td>&rho;<sub>resin</sub></td><td>1.26 g/cm<sup>3</sup></td></tr><tr><td>&rho;<sub>fiber</sub></td><td>1.80 g/cm<sup>3</sup></td></tr><tr><td>V<sub>f0</sub></td><td>57%</td></tr><tr><td>&lambda;<sub>11</sub></td><td>5.80 W/mK</td></tr><tr><td>&lambda;<sub>22</sub></td><td>0.57 W/mK</td></tr><tr><td>&lambda;<sub>33</sub></td><td>0.57 W/mK</td></tr></table>

Coefficients of Equation (7) and Equation (8) are defined in Table 1. Observe that the permeability depends on the fiber volume fraction, which in turn, depends on the pressure field, making the consolidation problem nonlinear itself. However, the permeability coefficients can be assumed to be uniform all over the domain, i.e. they do not depend on the space coordinates. This fact will allow for a very simple parametrization of the permeability tensor, see Section 3.3 for details.

## 3 Simulation based on reduced order modeling

In this section, we describe the simulation strategy, including the different simulation modules based on ROM techniques and their coupling. In particular, we recall the main features of the ROM methods that were implemented, emphasizing both the construction of the reduced models in the offline stage and their utilization in the online stage, as a part of the *Simulation App*.

The coupled model described in Section 2 involves three primary unknown fields, the temperature $T(x, t)$, the pressure $P(x, t)$ and the curing degree $\alpha(x, t)$. The viscosity field, $\eta(x, t)$, can be obtained as a post-processing of the temperature and curing degree fields. The nonlinearity associated to the coupling can be efficiently addressed by using a semi-implicit incremental time integration, and treating all coupling terms explicitly, therefore allowing for decoupling of the different problems. Thanks to this strategy, we are able to define different simulation modules, described below, one for each physics involved in the
manufacturing process. This approach is particularly effective in the context of ROM, because it allows applying the most suitable techniques depending on the nature of the equations.

It is to be noted that domain changes due too the applied compression will be neglected, i.e. the mould thickness reduction remains small enough.

## 3.1 Reduced order modeling methods

The simulation strategy combines three different numerical techniques: (i) the Proper Orthogonal Decomposition (POD), (ii) the Discrete Empirical Interpolation Method (DEIM) and (iii) the Proper Generalized Decomposition (PGD). In this section we summarize the main ingredients of these three techniques.

### 3.1.1 The Proper Orthogonal Decomposition

POD extracts the most significant components, measured in 2-norm, in the solutions of a parametric problem from the analysis of a set of “snapshots” (previously computed solutions of a given problem at different times and for different values of the model parameters) and uses them to approximate the solution for a new set of parameters up to a certain degree of accuracy.

Thus, POD allows expressing the unknown field involved in a generic problem $u(\mathbf{x}, t; \boldsymbol{\mu})$, with $\boldsymbol{\mu}$ the vector containing the model parameters, as

$$u(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_u} a_u^n(t; \boldsymbol{\mu}) \phi_u^n(\mathbf{x}), \qquad (9)$$

with usually much lower $N_u$ that the number of nodes, denoted by $I$, used for approximating the unknown field within the finite element framework. When considering discrete approximation spaces, the unknown vector at time $t_j$, $j = 1, \dots, J$, is denoted by $\mathbf{u}(t_j; \boldsymbol{\mu})$, that contains the value of the unknown field at each node $\mathbf{x}_i$, $i = 1, \dots, I$, can be expressed in the following reduced form:

$$\mathbf{u}(t_j; \boldsymbol{\mu}) = \mathbf{\Phi}_u \mathbf{a}_u(t_j; \boldsymbol{\mu}), \qquad (10)$$

defining the linear combination of vectors $\phi_u^n$, constituting the columns of matrix $\mathbf{\Phi}_u$, by the coefficients $\mathbf{a}_u(t_j; \boldsymbol{\mu})$, where the subscript $\cdot_u$ indicates that the reduced basis and the projection coefficients are associated to the (scalar) field $u$.

### 3.1.2 The Discrete Empirical Interpolation Method

When considering a generic nonlinear function $g(u(\mathbf{x}, t; \boldsymbol{\mu}))$ its evaluation in the reduced state space from the coefficients $\mathbf{a}_u(t_j; \boldsymbol{\mu})$ can be very expensive and, as discussed in [15, 43], it may compromise the efficiency of reduced model techniques. One possibility consists of using the POD just described to extract a reduced basis able to approximate the nonlinear term as

$$g(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_g} a_g^n(t; \boldsymbol{\mu}) \phi_g^n(\mathbf{x}). \qquad (11)$$
In order to compute the approximation coefficients $a_g^n$, with $n = 1, \dots, N_g$, the nonlinear function is reconstructed at any time $t_j$ at only $N_g$ nodes, denoted by $\mathbf{x}_g^i$, with $i = 1, \dots, N_g$. The subscript $\bullet_g$ stands for the fact of that points being related to the approximation of the nonlinear function on the reduced basis, $\phi_g^n$. It results in the following linear system to be solved each at each time instance and given set of parameters:

$$ \begin{pmatrix} g(\mathbf{x}_g^1, t_j; \boldsymbol{\mu}) \\ \vdots \\ g(\mathbf{x}_g^{N_g}, t_j; \boldsymbol{\mu}) \end{pmatrix} = \begin{pmatrix} \phi_g^1(\mathbf{x}_g^1) & \cdots & \phi_g^{N_g}(\mathbf{x}_g^1) \\ \vdots & \ddots & \vdots \\ \phi_g^1(\mathbf{x}_g^{N_g}) & \cdots & \phi_g^{N_g}(\mathbf{x}_g^{N_g}) \end{pmatrix} \begin{pmatrix} a_g^1(t_j; \boldsymbol{\mu}) \\ \vdots \\ a_g^{N_g}(t_j; \boldsymbol{\mu}) \end{pmatrix} .(12) $$

The remaining question concerns the choice of points $\mathbf{x}_g^i$, $i = 1, \dots, N_g$. These points are a priori arbitrary, the only constraint being that the matrix involved in Equation (12) be invertible. However, in [43] authors propose to define these points using a greedy strategy in order to locate them to capture as much information as possible. The resulting points were called “magic points” and for the sake of completeness we describe below their calculation. We start by considering

$$ \mathbf{x}_g^1 = \operatorname{argmax}_{\mathbf{x}} |\phi_g^1(\mathbf{x})|. \qquad (13) $$

Then we compute $d_1$ from

$$ d_1 \phi_g^1(\mathbf{x}_g^1) = \phi_g^2(\mathbf{x}_g^1), \qquad (14) $$

that allows defining the residual $r_2(\mathbf{x})$, i.e. a contribution in $\phi_g^2(\mathbf{x})$ that cannot be explained by $\phi_g^1(\mathbf{x})$, from which computing point $\mathbf{x}_g^2$:

$$ \mathbf{x}_g^2 = \operatorname{argmax}_{\mathbf{x}} |r_2(\mathbf{x})| \quad \text{with} \quad r_2(\mathbf{x}) = \phi_g^2(\mathbf{x}) - d_1 \phi_g^1(\mathbf{x}). \qquad (15) $$

As by construction $r_2(\mathbf{x}_g^1) = 0$ we can ensure $\mathbf{x}_g^2 \neq \mathbf{x}_g^1$. The procedure is generalized for obtaining the other points involved in the interpolation procedure. Thus, for obtaining point $\mathbf{x}_g^i$ we consider

$$ \mathbf{x}_g^i = \operatorname{argmax}_{\mathbf{x}} |r_i(\mathbf{x})| \quad \text{with} \quad r_i(\mathbf{x}) = \phi_g^i(\mathbf{x}) - \sum_{j=1}^{i-1} d_j \phi_g^j(\mathbf{x}). \qquad (16) $$

The coefficients $d_1, \dots, d_{i-1}$ must be chosen for ensuring that $\mathbf{x}_g^i \neq \mathbf{x}_g^j, \forall j < i$. For this purpose we enforce that the residual $r_i(\mathbf{x})$ vanishes at each location $\mathbf{x}_g^j$, with $j < i$ by solving:

$$ r_i(\mathbf{x}_g^j) = 0 = \phi_g^i(\mathbf{x}_g^j) - \sum_{l=1}^{i-1} d_l \phi_g^l(\mathbf{x}_g^l), \quad j = 1, \dots, i-1, \qquad (17) $$

that constitutes a linear system whose solution leads to the sought-after coefficients $d_1, \dots, d_{i-1}$.
### 3.1.3 The Proper Generalized Decomposition

Most of the existing model reduction techniques proceed by extracting a suitable reduced basis and then projecting the problem solution on it. Thus, the reduced basis construction precedes its use in the solution procedure, and one must be careful on the suitability of a particular reduced basis when employed for representing the solution of a particular problem.

This issue disappears if the approximation basis is constructed at the same time that the problem is solved. Thus, each problem has its associated basis in which its solution is expressed. One could consider few terms in its approximation, leading to a reduced representation, or all the terms needed for approximating the solution up to a certain accuracy level. The Proper Generalized Decomposition (PGD) proceeds in this manner.

When addressing the solution of a parametric problem $u(x, t; \mu)$ its PGD-based separated representation, including parameters as extra-coordinates, reads

$$u(x, t, \mu) \approx \sum_{n=1}^{N_u} \phi_x^n(x) \phi_t^n(t) \phi_\mu^n(\mu). \qquad (18)$$

For additional details the interested reader can refer to [38] and the numerous references therein.

## 3.2 The thermo-kinetic simulation module

A POD approach is applied in order to reduce the computational complexity of the thermo-kinetic simulation. Reduced basis are computed for both primary fields, temperature and curing degree. Other reduced basis are also computed in order to approximate nonlinearities using the DEIM. The training stage, i.e. snapshots generation, and the reduced basis extraction are explained in Section 3.2.1. Then, in Section 3.2.2, the ROM for the thermo-kinetic simulation is assembled. The reduced version of the heat equation is formed by standard Galerkin projection onto the temperature reduced basis. The kinetic equation can be treated in a different manner thanks to its local nature. The computational complexity is reduced in this case by integrating Equation (1b) only in a well-chosen subset of mesh nodes; in particular, that subset will be computed using DEIM.

### 3.2.1 Training stage and reduced basis extraction

In the training stage, a series of simulations were carried out with a full-order solver based on standard numerical techniques, for different choices of the model parameters. In our case, the following model parameters were identified:

* Initial degree of curing, $\alpha_0$, that may differ from one OGV part to another depending on the time the pre-pregs stay out of the freezer before manufacturing.
* Initial temperature, $T_0$, that may change depending on the ambient conditions.

* Temperature cycle imposed at both top and bottom walls of the mold, denoted by $T_t(t)$ and $T_b(t)$, respectively. In the actual process, however, a fixed temperature is kept all along the process, i.e. the variability at both the beginning and the end can be neglected. In addition, both top and bottom temperatures are the same in practice so we can simply denote $T_c \equiv T_t = T_b$.

* Initial fiber volume fraction, $V_{f0}$, that may change if the ply stack is modified by adding or removing plies.

We denote by $\boldsymbol{\mu}_m = (\alpha_0, T_0, T_c, V_{f0})_m$, with $1 \le m \le M$, the parameter sampling in order to compute the snapshots. Recall that the objective is to train a ROM from these snapshots. The nonlinear thermo-kinetic coupled problem given by Equation (1) was therefore solved to obtain the space-time temperature and curing degree space-time fields, i.e. $T(\mathbf{x}, t; \boldsymbol{\mu}_m)$ and $\alpha(\mathbf{x}, t; \boldsymbol{\mu}_m)$.

From these collected data, a reduced basis is extracted for both temperature and curing degree, allowing the introduction of the following approximations:

$$T(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_T} a_T^n(t; \boldsymbol{\mu}) \phi_T^n(\mathbf{x}) \quad \text{and} \quad \alpha(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_\alpha} a_\alpha^n(t; \boldsymbol{\mu}) \phi_\alpha^n(\mathbf{x}), \quad (19)$$

where $\phi_T^n$ and $\phi_\alpha^n$ denote the temperature and curing degree basis functions, respectively. Coefficients $a_T^n$ are to be determined by solving a reduced system, which will be formed by Galerkin projection onto the reduced basis $\phi_T^n$. On the other hand, coefficients $a_\alpha^n$ are to be determined by integrating the kinetic equation only in a subset of mesh nodes, denoted by $\{\mathbf{x}_\alpha^n\}_{n=1}^{N_\alpha}$. This means that the kinetic equation may be integrated at only $N_\alpha$ locations (very few in practice) leading to important computational time savings. With the reduced basis for the curing degree at hand, these nodes are chosen straightforwardly as explained in Section 3.1.2.

The nonlinear terms in the thermo-kinetic model, namely the specific heat and the curing function, also need a special treatment in order to build an efficient reduced solver. According to Section 3.1.2, a reduced basis has to be formed for both specific heat and curing function.

Observe that the reduced basis for the curing degree can also be used to approximate the curing function. This is because the snapshots of the curing rate, $\dot{\alpha}(\mathbf{x}, t; \boldsymbol{\mu}_m)$ (equivalently, $f^{\text{kin}}(\mathbf{x}, t; \boldsymbol{\mu}_m)$), are linear combination of the snapshots of the curing degree $\alpha(\mathbf{x}, t; \boldsymbol{\mu}_m)$. Therefore, it is clear that the space spanned by both basis would be the same, and in consequence the curing function is approximated using the curing degree basis:

$$f^{\text{kin}}(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_f} a_f^n(t; \boldsymbol{\mu}) \phi_f^n(\mathbf{x}) \quad \text{with} \quad \phi_f^n(\mathbf{x}) \equiv \phi_\alpha^n(\mathbf{x}) \quad (N_f = N_\alpha). \quad (20)$$
As a corollary, the DEIM points for approximating the curing function will also be the same than those for the curing degree, i.e. {$\mathbf{x}_f^n$}$_{n=1}^{N_f}$ ≡ {$\mathbf{x}_\alpha^n$}$_{n=1}^{N_\alpha}$. Coefficients $a_f^n$ are to be determined by solving an interpolation problem, as shown in Section 3.1.2.

Regarding the specific heat approximation, the snapshots $C_p(\mathbf{x}, t; \boldsymbol{\mu}_m)$ are obtained simply by evaluating Equation (4) using the temperature and curing degree snapshots. Therefore, we introduce the following approximation:

$$C_p(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_C} a_C^n(t; \boldsymbol{\mu}) \phi_C^n(\mathbf{x}), \qquad (21)$$

where $\phi_C^n$ denote the specific heat basis functions. With the reduced basis at hand, the corresponding DEIM points can be computed; we denote them by {$\mathbf{x}_C^n$}$_{n=1}^{N_C}$. Coefficients $a_C^n$ are to be determined by solving an interpolation problem, as explained in Section 3.1.2.

### 3.2.2 Assembling the Reduced order model

Let us first focus on the heat equation. We first introduce the approximations of the nonlinear terms, i.e. the specific heat and kinetic rate, given by Equation (20) and Equation (21):

$$\rho \left\{ \sum_{n=1}^{N_C} a_C^n(t; \boldsymbol{\mu}) \phi_C^n(\mathbf{x}) \right\} \frac{\partial T}{\partial t} = \lambda \nabla^2 T + \rho \Delta H \left\{ \sum_{n=1}^{N_\alpha} a_f^n(t; \boldsymbol{\mu}) \phi_\alpha^n(\mathbf{x}) \right\}. \quad (22)$$

**Remark 2** (On the approximation of the kinetic rate). As explained in Section 3.2.1, the kinetic rate is approximated using the basis for the curing degree, and so in the right-hand side of Equation (22) we use $\phi_\alpha^n$. However, it is worth to remark for the sake of clarity that we keep the coefficients $a_f^n$, which are obviously not the same than $a_\alpha^n$.

Using a standard finite element discretization on the mesh, it results in the following discrete matrix form:

$$\left\{ \sum_{n=1}^{N_C} a_C^n(t; \boldsymbol{\mu}) \mathbf{M}_n \right\} \dot{\mathbf{T}} + \mathbf{K} \mathbf{T} = \left\{ \sum_{n=1}^{N_\alpha} a_f^n(t; \boldsymbol{\mu}) \mathbf{b}_n \right\}, \quad (23)$$

where $\mathbf{M}_n$ are the different mass matrices, each one related to a specific heat basis function, $\phi_C^n(\mathbf{x})$; $\mathbf{K}$ is the conductivity matrix (not to be confused with the permeability tensor, although no distinct notation is used); and $\mathbf{b}_n$ are the different source vectors related to each curing function basis $\phi_\alpha^n$.

By performing the Galerkin projection of Equation (23) onto the reduced basis for the temperature, denoted in its discrete form by $\mathbf{\Phi}_T$, and taking into account that discrete equivalent of Equation (19) is $\mathbf{T} = \mathbf{\Phi}_T \mathbf{a}_T$, we get:

$$\left\{ \sum_{n=1}^{N_C} a_C^n(t; \boldsymbol{\mu}) \tilde{\mathbf{M}}_n \right\} \dot{\mathbf{a}}_T + \tilde{\mathbf{K}} \mathbf{a}_T = \left\{ \sum_{n=1}^{N_\alpha} a_f^n(t; \boldsymbol{\mu}) \tilde{\mathbf{b}}_n \right\}, \quad (24)$$
with the following definitions for the reduced operators: $\tilde{M}_n = \Phi_T^t M_n \Phi_T \in \mathbb{R}^{N_T \times N_T}$, $\tilde{K} = \Phi_T^t K \Phi_T \in \mathbb{R}^{N_T \times N_T}$ and $\tilde{b}_n = \Phi_T^t b_n \in \mathbb{R}^{N_T \times 1}$. The transpose is denoted by $\bullet^t$.

Recall that coefficients $a_C^n$ and $a_f^n$ are computed, at each time, by solving a linear system of size $N_C$ and $N_\alpha$, respectively, as explained in Section 3.1.2.

Finally, the reduced version of the kinetic equation is trivial thanks to its locality. In order to compute the coefficients $a_\alpha^n$ at each time, we are integrating Equation (1b) at the DEIM points related to the curing degree:

$$ \dot{\alpha}(\mathbf{x}_\alpha^n, t; \boldsymbol{\mu}) = f^{\text{kin}}(T(\mathbf{x}_\alpha^n, t; \boldsymbol{\mu}), \alpha(\mathbf{x}_\alpha^n, t; \boldsymbol{\mu})) \quad \text{with} \quad 1 \le n \le N_\alpha, \quad (25) $$

which of course gives $N_\alpha$ uncoupled ODEs. In discrete form:

$$ \dot{\tilde{\alpha}}(t; \boldsymbol{\mu}) = f^{\text{kin}}(\tilde{T}(t; \boldsymbol{\mu}), \tilde{\alpha}(t; \boldsymbol{\mu})), \quad (26) $$

where $\dot{\tilde{\alpha}}, \tilde{T}, \tilde{\alpha} \in \mathbb{R}^{N_\alpha \times 1}$. Coefficients $a_\alpha^n(t; \boldsymbol{\mu})$ in Equation (19) are computed by solving the following interpolation problem:

$$ \alpha(\mathbf{x}_\alpha^m, t; \boldsymbol{\mu}) = \sum_{n=1}^{N_\alpha} a_\alpha^n(t; \boldsymbol{\mu}) \phi_\alpha^n(\mathbf{x}_\alpha^m) \quad \text{with} \quad 1 \le m \le N_\alpha. \quad (27) $$

### 3.2.3 Implementation details and results

We report hereafter the details for the numerical implementation of the full-order solver, used for computing the snapshots:

* The Finite Element mesh is composed of 116, 136 HEX8 elements, with 147, 245 nodes.

* A first-order semi-implicit integration scheme was applied. The simulated time is 2 hours of process, with 20, 000 time steps.

* A Preconditioned Conjugate Gradient is used as linear solver.

* The typical running time, for a single parameter execution, is 5 hours on a 64bit machine with the following specifications: 2.9 GHz Intel Core i5 with 16 Gb 1867 MHz DDR3, running under MacOS X 10.10.5.

Concerning the snapshots generation, we considered random perturbations around the nominal parameters: $\mu_{\text{nominal}} = (\alpha_0 = 0, T_0 = 288 \text{ K}, T_c = 438 \text{ K}, V_{f0} = 57\%)$. The amplitude of the perturbations is in the following ranges: $\alpha_0 \in [0, 0.05]$ (+5%), $T_0 \in [273, 303] \text{ K}$ ($\pm 5\%$), $T_c \in [416, 460] \text{ K}$ ($\pm 5\%$), $V_{f0} \in [50, 64]\%$ ($\pm 7\%$). Reduced basis of the following dimension (i.e. number of degrees of freedom of the ROM) were found: $N_T = 10$, $N_\alpha = 7$ and $N_C = 5$, which proves the pertinency of using low dimensional approximations for the thermo-kinetic model. We report hereafter the details for the numerical implementation of the ROM solver:

* A variable order (1-5) Numerical Differentiation Formula (MATLAB ode15s) time integrator was used. The simulated time is 2 hours and the number of time steps varies adaptively.
<figure><img src="image_6.png" /><figcaption>Figure 4: First four most significant temperature modes: $\phi_T^1(\mathbf{x})$ (top-left); $\phi_T^2(\mathbf{x})$ (top-right), $\phi_T^3(\mathbf{x})$ (bottom-left) and $\phi_T^4(\mathbf{x})$ (bottom-right).</figcaption></figure>

* The typical running time is in the range 1 – 10 seconds depending on the process conditions, using the same machine described previously. It is important to emphasize that the simulation took few seconds to completion instead of the many hours needed when solving the full order problem.

* The random access memory (RAM) consumption was measured by running 5 simulations, each one of them for a different set of process parameters selected randomly. The simulation time was set to 1 hour. In such conditions, the averaged memory required is 236 Mb. More details about the RAM consumption of the final implementation will be given in Section 4.

* A posteriori validation of the reduced model was carried out by comparing results with full-order FEM simulations for sets of randomly generated process parameters. The observed relative error measured in the maximum norm is typically of the order of $10^{-5}$.

Figure 4 and Figure 5 depict the first four modes (the most significant ones) of the temperature and the curing degree. Figure 6 shows the time evolution of the weights related to the temperature and degree of curing modes.

Figure 7 represents the temperature solution issued from the reduced simulation, at four different times, after recombining the basis functions with the reduced coordinates.
<figure><img src="image_7.png" /><figcaption>Figure 5: First four most significant curing degree modes: $\phi_{\alpha}^{1}(\mathbf{x})$ (top-left); $\phi_{\alpha}^{2}(\mathbf{x})$ (top-right), $\phi_{\alpha}^{3}(\mathbf{x})$ (bottom-left) and $\phi_{\alpha}^{4}(\mathbf{x})$ (bottom-right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Time evolution of the temperature (left) and curing degree (right) reduced coordinates, $a_{T}^{n}(t; \boldsymbol{\mu})$, related to the first five temperature modes (i.e. $n \le 5$) during the first 100 sec of the simulation, corresponding to a randomly chosen parameter set.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Temperature fields at four different times from the thermo-kinetic ROM: $t = T/4$ (top-left); $t = T/2$ (top-right); $t = 3T/3$ (bottom-left) and $t = T$ (bottom-right).</figcaption></figure>

## 3.3 The consolidation simulation module

The PGD method was applied in order to reduce the computational complexity of the consolidation model. Recall that the consolidation model is nonlinear itself, because the permeability tensor depends on the fiber volume fraction. However, as explained in Section 2.3, the permeability tensor can be assumed uniform, i.e. not varying through the domain. In consequence, observe that the problem could be linearized by simply considering the permeability tensor as an extra-coordinate in the PGD framework. A similar approach was already addressed in the context of nonlinear soil dynamics [44]. Hence, in this scenario, we would have a parametric pressure field in the form: $P(x, K)$. And more importantly, the computation of such parametric would now involve a linear problem. This can be understood by thinking of the parametric solution as being a linear solver inside a fixed-point nonlinear scheme.

However, there is an additional difficulty concerning the parametrization of the viscosity field. A low-dimensional parametrization of such field is sought by following a similar approach to that one explained in Section 3.2. In fact, we are seeking for an approximation of the inverse of the viscosity field, see Equation (5), as follows:

$$ \eta^{-1}(\mathbf{x}, t; \boldsymbol{\mu}) \approx \sum_{n=1}^{N_{\eta}} a_{\eta}^{n}(t; \boldsymbol{\mu}) \phi_{\eta}^{n}(\mathbf{x}), \qquad (28) $$

where $\phi_{\eta}^{n}$ denotes the basis functions for the inverse of the viscosity, computed from snapshots $\eta^{-1}(\mathbf{x}, t; \boldsymbol{\mu}_{m})$, $1 \le m \le M$, by simple evaluation of Equation (7).
Observe that coefficients $a_{\eta}^n$ can be seen in fact as the coordinates of the inverse of the viscosity on the reduced basis, at every instant and choice of parameters. We are therefore introducing them as extra-coordinates in the PGD parametric solution, i.e. $P(\mathbf{x}, \mathbf{K}, \mathbf{a}_{\eta})$. Such parametric solution is now able to provide the pressure field, not only for every possible permeability tensor, but also for every viscosity field living in the reduced basis $\phi_{\eta}^n$. Of course, we hope $N_{\eta}$ to be small enough.

We shall explain in next section how to compute the coefficients $\mathbf{a}_{\eta}$ in order to be able to evaluate the parametric solution.

**Remark 3** (On the dependence of the pressure field on the closure rate). *Simulating the closure cycle of the press requires obtaining the pressure field, for a given viscosity field, permeability tensor (i.e. fiber volume fraction) and an imposed closure rate. However, it is worth to remark that the closure rate does not need to be considered as an extra-coordinate of the parametric solution, because of the linearity of the problem to be solved. In other words, the pressure field depends linearly on the closure rate.*

The parametric solution of the flow problem, Equation (29), within the PGD framework to obtain a separated representation of the parametric pressure field:

$$P(\mathbf{x}, \mathbf{K}, \mathbf{a}_{\eta}) = \sum_{n=1}^{N_P} \phi_x^n(\mathbf{x}) \phi_K^n(\mathbf{K}) \phi_a^n(\mathbf{a}_{\eta}). \quad (29)$$

Even if such solution is related to a multidimensional problem, its linearity makes possible to compute it very efficiently by means of the PGD.

### 3.3.1 Implementation details and results

We report hereafter the details corresponding to the calculation of the PGD parametric pressure solution:

* The same mesh already described in Section 3.2.3 was used.

* The inverse of the viscosity is parametrized using $N_{\eta} = 3$ coordinates, corresponding to the first three modes $\phi_{\eta}^1(\mathbf{x})$, $\phi_{\eta}^2(\mathbf{x})$ and $\phi_{\eta}^3(\mathbf{x})$ that can represent the 95% of the solution.

* The FE mesh for the viscosity in the phase-space $(a_{\eta}^1, a_{\eta}^2, a_{\eta}^3)$ is made of 6,307 TET4 elements and 1,387 nodes.

* The FE mesh for the permeability phase-space $k_{11}, k_{22}, k_{33}$ is made of EDGE2 elements $1D \times 1D \times 1D$, 100 nodes each.

* The PGD convergence criterion was set to $10^{-3}$ on the relative norm of the residual.

Typical execution times are not reported for the consolidation ROM because they are completely negligible: evaluating the parametric pressure is straightforward.
<figure><img src="image_10.png" /><figcaption>Figure 8: First four most significant inverse viscosity modes: $\phi_{\eta}^{1}(\mathbf{x})$ (top-left); $\phi_{\eta}^{2}(\mathbf{x})$ (top-right), $\phi_{\eta}^{3}(\mathbf{x})$ (bottom-left) and $\phi_{\eta}^{4}(\mathbf{x})$ (bottom-right).</figcaption></figure>

Figure 8 depicts the first four modes (the most significant ones) of the inverse of the viscosity respectively, although only three were kept as explained before. The time evolution of the three associated coefficients, $a_{\eta}^{1}(t)$, $a_{\eta}^{2}(t)$ and $a_{\eta}^{3}(t)$, describes a closed curve in the 3D espace defined by them, as depicted in Figure 9, since both the initial and final states are characterized by a uniform high viscosity field.

Finally Figure 10 shows for four different rheologies expressed from vector $\mathbf{a}_{\eta}$ in the viscosity phase-space as well as their four associated pressure fields.

## 3.4 Online simulation strategy: simulation modules coupling

We summarize in this section the operation of the reduced models, with the help of Figure 11. Suppose that the reduced coordinates of the temperature and the curing degree are known from the previous time step; we denote $\mathbf{a}_{T}(t_{j-1})$ and $\mathbf{a}_{\alpha}(t_{j-1})$, respectively. Note that at time $t_0$, these reduced coordinates are obtained from the initial conditions for both temperature and curing degree. Then:

1. We evaluate both the specific heat and the curing function at the previous time step, i.e. a semi-implicit scheme is used. The reduced coordinates of both specific heat and curing function, denoted by $\mathbf{a}_{C}(t_{j-1})$ and $\mathbf{a}_{f}(t_{j-1})$ respectively, are computed using DEIM.

2. With those reduced coordinates at hand, the thermo-kinetic ROM, Equa-
<figure><img src="image_11.png" /><figcaption>Figure 9: Trajectory of coefficients $a_{\eta}^n(t)$ (left) and parametric domain ensuring positive viscosities (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Four examples showing arbitrary viscosities in the phase-space and their corresponding pressure field. For each example: (i) pressure field in the OGV part (left); (ii) reduced coordinates of the viscosity (red point) in the positive viscosities phase space (right).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Schematics showing the operation of both thermo-kinetic and consolidation simulation modules based on ROM, as well as their coupling.</figcaption></figure>
tion (24) and Equation (26), can be assembled as explained in Section 3.2.

3. By performing a time increment, the reduced coordinates of both temperature and curing degree can be obtained at time $t_j$.

4. Then, the reduced coordinates of the viscosity, $a_η$, can be computed by evaluating the chemo-rheology model and using DEIM.

5. With the previous information at hand, the consolidation module can be run. It involves solving a nonlinear problem because of the press control system, see Section 2.1. The nonlinear iterations are indexed by $ℓ$. We start by assuming that the closure rate is the same than in the previous time increment, i.e. $\dot{U}^{\ell=0}(t_j) = \dot{U}(t_{j-1})$. Similarly, we consider $V_f^{\ell=0}(t_j) = V_f(t_{j-1})$.

6. Using Equation (8), the permeability tensor can be computed.

7. With the permeability tensor and the reduced coordinates of the viscosity, the parametric PGD solution can be particularized. Recall that this solution has been computed for unitary velocity and so it must be scaled by the actual velocity. We denote by $\mathbf{P}^\ell(t_j)$ the pressure field at time $t_j$ and iteration $ℓ$.

8. The velocity field can be then computed using Equation (5). Integrating on the free boundary $S$, see Figure 3, the volume of resin drained can be computed, $\Delta V_r^\ell$. Note that knowing the loss of resin it is trivial to compute the current fiber volume fraction.

9. On the other hand, from the pressure field it is possible to obtain the vertical component of the reaction force to be applied by the press so as to maintain the actual closure rate. The press force can be computed by integrating the pressure field on the upper boundary $S^+$, see Figure 3. However, note that this force only represents the fluid part, i.e. the resin, but it does not take into account the fiber contribution. In order to account for this extra contribution, we consider the Kim's model.

10. With the press force, $F^\ell(t_j)$ at hand, the press control can be checked. A Newton-Raphson algorithm is used in order to compute both closure rate and fiber volume fraction correction, denoted by $\Delta \dot{U}$ and $\Delta V_f$.

11. From the corrections, both closure rate and fiber volume fraction can be updated. At convergence, we set: $\dot{U}(t_j) \leftarrow \dot{U}^{\ell*}(t_j)$ and $V_f(t_j) \leftarrow V_f^{\ell*}(t_j)$, assuming the algorithm converges after $\ell^*$ iterations.

12. Observe that the semi-implicit linearization allows excluding the thermo-kinetic module from the nonlinear problem. Otherwise, at each nonlinear iteration, it would be required to come back to the thermo-kinetic module, to recompute the specific heat, to solve the thermo-kinetic ROM, and so on.
## 4 A Simulation App for the OGV manufacturing process

The *Simulation App* for the OGV manufacturing process is a process-specific application that allows the user to simulate almost in real time different process conditions and visualize the simulation results. Here we understand by real time simulation the one able to provide the results with no perceptible delay after the user makes its request, i.e. in the order of half a second. However, since the *Simulation App* is designed to be used directly in the manufacturing facility, the input and output delays are also relevant. As it will be explained later, both input and output interfaces were designed so as to enhance the user's reactivity by restricting the input data to the bare minimum, while only relevant output information is displayed. Typically, a standard user takes no more than one minute in order to set up a new simulation, while the time consumed in visualizing and interpreting the results depends mostly on the user's knowledge and experience.

The user interacts with the *Simulation App* through a basic Graphics User Interface (GUI). In order to demonstrate the feasibility and potentiality of this kind of applications, we developed a demonstration version in the MATLAB® environment using GUIDE under MacOS X 10.10.5. This allows creating the GUI very easily. The MATLAB Application Compiler™ was used in order to create a standalone executable file that could run outside MATLAB® and be easily transferred to the final user for demonstration purposes. The standalone version of the *Simulation App* has been tested on Microsoft Windows® 7 OS. Other proper implementations are of course possible although not explored in this paper, as they are outside of its scope.

The concept of a *Simulation App* offers several potentialities. Since the application is process-specific, rather than implementing a general purpose visualization environment, it is possible to first identify the set of quantities of interest as output and then implement a simple and specific visualization interface. For example, if we know in advance that the maximum pressure gradient is an important indicator for defects, as it is the case in OGV manufacturing processes, we can include a functionality that displays the maximum pressure gradient and its location, at each time step. Then, the user can access to this information by simply activating a checkbox.

A similar discussion can be addressed regarding the data input. In a general purpose simulation software, a number of simulation parameters must be entered as user input. Although these have a significant impact within the modeling and simulation framework, they have little significance for the process designer. In order to get to these, the user must convert process parameters (physical measures, material properties) into simulation parameters (boundary and initial conditions, model adjustments) before running the simulation. The *Simulation App* is designed to perform this conversion automatically according to rules provided beforehand. For instance, the initial fiber volume fraction is computed from process parameters, such as the following: area baseline plykit
<figure><img src="image_14.png" /><figcaption>Figure 12: *Simulation App* pre-processing tabs: Parameters tab (left), Data Tab (center) and Simulation Tab (right).</figcaption></figure>

(stack of composite layers), the resin and fiber mass density, the mass of the convex and concave stack and the fiber areal weight (FAW) and the volatile content, see Figure 12. Therefore, the process designer only enters process parameters, which are the real meaningful information to them, without the need of expertise in multi-physics modeling or numerical methods. The difference between process and simulation parameters is kept transparent to the user in the GUI. By consequence, a *Simulation App* is not only simulation tool but it also integrates the process knowledge via process-specific inputs and outputs.

The *Simulation App* is composed of three modules:

1. *Pre-processing*. This module performs two basic operations: data loading and parameters conversion. Data loading simply reads all pre-computed data required to run the reduced model. This operation only needs to be done once after launching the application, and it is performed in the GUI via the Load Data button. Parameters conversion gathers data entered by the user and computes the simulation parameters. This operation is performed in the GUI via two different tabs: the Parameters Tab, in which some default values that normally do not change are proposed (e.g. resin and fiber specific weight), and the Data Tab, in which process parameters are defined (e.g. temperature cycle, closure rate, etc.) (see Figure 12). The conversion operation is launched by the Update button, and the initial degree of curing and the initial fiber volume content, both being simulation parameters, are computed.
2. *Simulation*. This module is driven by a principal function that governs the interaction between the two reduced models explained in previous sections: POD reduced model for the thermo-kinetic simulation and PGD reduced model for the consolidation simulation. Basically, this module takes all data defined in the pre-processing module and runs the ROM model in the time interval of interest defined by the user, as explained in Section 3.4. The user can also choose the number of equally spaced time frames at which the solution wants to be accessed. Additionally, the user can also demand to access to the solution at particular time frames of interest, such as minimum viscosity time frame, see Figure 12 for details.

3. *Visualization*. This module allows accessing to both field data and quantities of interest, at any desired time frame. This operation is performed by the GUI via two different tabs: the Display Tab and the Quantities of Interest Tab. In the Display Tab, seven scalar fields can be visualized: curing degree, temperature, viscosity, pressure and the three components of the pressure gradient. They are simultaneously displayed on the external boundary of the part as well as in five sections, see Figure 13. The section view is necessary to appreciate, for instance, the temperature evolution inside the part, which may be higher than in the external boundary due to the exothermic kinetics. Additionally, the location and magnitude of the maximum and minimum of each field can be visualized. In the Quantities of Interest Tab, the time evolution of nine process indicators can be visualized: closure rate, press force, fiber volume fraction or part mass, for instance. See Figure 14 for details. These are the real meaningful information upon which a process designer can evaluate if the process setting is operating as desired or not.

The RAM consumption of the MATLAB® implementation of the *Simulation App* was measured. Opening the application requires 239 Mb, whereas the amount of memory increases up to 504 Mb after loading data. The amount of memory required for running the simulation depends on the simulation timespan as well as on the process parameters. In the conditions already described in Section 3.2.3, 236 Mb are required on average, which brings the total to 740 Mb. However, the most demanding operation in terms of memory consumption is the Visualization module, which requires up to 1160 Mb. In general, the memory consumption could be drastically reduced in a proper (non demonstrator) implementation of the application. Some parts of the code could also be optimized, for instance, by avoiding reconstruction of the entire fields, since only some slices and the external surface are visualized.

ROM-based *Simulation Apps* can be a powerful tool for complex composite processes to increase the entitlement yield by adapting for the variation that comes from material chemistry and physical properties in addition to thermal and pressure histories applied during the process. Typical entitlement yield is limited because of the inherent variations and multi-physics interactions. Part quality loss is due to either (i) internal defects, (ii) not meeting dimensional
<figure><img src="image_15.png" /><figcaption>Figure 13: *Simulation App* post-processing: Display tab. Curing degree field after 3.05 min (top). Temperature field at after 20.30 min (mid). Pressure field after 60 min (bottom).</figcaption></figure>
requirements or (iii) poor internal fiber matrix structure. A successful man-
ufacturing process must minimize all three quality components. The *Simulation App* can be implemented in the manufacturing process seamlessly to make
real-time decisions regarding process adaptation possible overcoming inherent
variation and truly increasing the entitlement yield of the process. As an ex-
ample, the *Simulation App* will take real process measurements as inputs which
capture the incoming variation in the pre-processing section and provide quan-
tities of interest e.g. maximum pressure gradient that can be relate to quality.
If the quality does not meet requirements, pre-computed sensitivity results can
be used to identify corrective process change to bring it within requirements.
The corrective process change can be as simple as choosing among predefined
process cycles. All this can be made possible only because of real-time models
capable of reflecting physics of the entire process.

5 Conclusions

In this paper we have introduced the *Simulation App* concept, a process-specific simulation tool based on reduced order modeling techniques. Using a combination of them, we were able to reduce from several hours to few seconds the computational time of a coupled multi-physics and strongly nonlinear model describing the manufacturing process of a composite outlet guide vane. Defining a coupling strategy between the different reduced order models was an essential part of this work. In particular, we presented an approach based on an appropriate parametrization of the coupling fields. The use of such fast simulation in a real-time decision making environment being possible, process-specific pre-processing and visualization functionalities were added, leading to the *Simulation App* concept.

It has also been demonstrated that the *Simulation App* provides several
advantages over general-purpose simulation software, especially if simulation
wants to be used directly in the manufacturing facility. In addition to the com-
putational time reduction, the process specificity of the *Simulation App* makes
it possible to conceive simple yet functional graphic interfaces, for both data
input and visualization. The process designer is asked to enter only process
parameters while the simulation parameters are limited to the bare minimum.
Similarly, the visualization module was designed so as to display only the rele-
vant information, mainly the process indicators upon which decisions are made.
Note that, by definition, process indicators are specific to a particular pro-
cess, which is contrary to the spirit of general-purpose software. Therefore, the
*Simulation App* establishes a link between process parameters and process in-
dicators through a comprehensive numerical simulation which includes not only
the physics and their couplings but also technological constraints such as the
control loop of the process, if any.

Finally, even if the *Simulation App* has been presented on a case study of
industrial interest, it is easy to imagine building similar tools for other manu-
facturing processes involving different physics, materials and technology.
<figure><img src="image_16.png" /><figcaption>Figure 14: *Simulation App* post-processing: Quantities of Interest tab. Closure rate (top), fiber volume rate (mid) and part mass evolution (bottom) during the simulated time, 60 min.</figcaption></figure>
References

[1] S.R. White, H.T. Hahn. Cure cycle optimization for the reduction of processing-induced residual stresses in composite materials. Journal of Composite Materials, 27(14), 1352-1378, 1993.

[2] V.K. Pillai, A.N. Beris, P. Dhurjati. Intelligent curing of thick section composites using a knowledge-based system. Journal of Composite Materials, 31(1), 2251, 1997.

[3] A. Shojaei, S.R. Ghaffarian, S.M.H. Karimian. Modeling and simulation approaches in the resin transfer molding process: a review. Polymer Composites, 24(4), 525-544, 2003.

[4] C.H. Park, L. Woo. Modeling void formation and unsaturated flow in liquid composite molding processes: a survey and review. Journal of reinforced plastics and composites, 30(11), 957-977, 2011.

[5] J.L Kardos, M.P. Dudukovic, R. Dave. Void growth and resin transport during processing of thermosettingmatrix composites. In Epoxy resins and composites IV, 101-123, Springer Berlin Heidelberg, 1986.

[6] N. Zobeiry, R. Vaziri, A. Poursartip. Computationally efficient pseudo-viscoelastic models for evaluation of residual stresses in thermoset polymer composites during cure. Composites Part A: Applied Science and Manufacturing, 41(2), 247-256, 2010.

[7] A.L. Mounier, C. Binetruy, P. Krawczak. Multipurpose carbon fiber sensor design for analysis and monitoring of the resin transfer molding of polymer composites. Polymer Composites, 26(5), 717730, 2005.

[8] E. Schmachtenberg, J. zur Heide Schulte, J. Topker. Application of ultrasonics for the process control of resin transfer moulding (RTM). Polymer Testing, 24(3), 330338, 2005.

[9] J.M. Lawrence, K.T. Hsiao, R.C. Don, P. Simacek, G. Estrada, E.M. Sozer, H.C. Stadtfeld, S.G. Advani. An approach to couple mold design and online control to manufacture complex composite parts by resin transfer molding. Composites Part A: Applied Science and Manufacturing, 33(7), 981990, 2002.

[10] F. Chinesta, R. Keunings, A. Leygue. The Proper Generalized Decomposition for advanced numerical simulations. A primer, Springerbriefs, Springer, 2014.

[11] A. Antoulas, D. Sorensen, S. Gugercin. A survey of model reduction methods for large-scale systems. Contemp. Math., 280, 193220, 2001.

[12] A. Quarteroni, A. Manzoni, F. Negri. Reduced Basis Methods for Partial Differential Equations: An Introduction. Springer International Publishing, 2015.
[13] S. Volkwein. Model reduction using proper orthogonal decomposition. Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz, 2011, http://www.uni-graz.at/imawww/volkwein/POD.pdf

[14] K. Willcox, J. Peraire. Balanced model reduction via the proper orthogonal decomposition. AIAA journal, 40(11), 2323-2330, 2002.

[15] S. Chaturantabut, D.C. Sorensen. Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput., 32, 2737-2764, 2010.

[16] B. Peherstorfer, D. Butnaru, K. Willcox, H.J. Bungartz (2014). Localized discrete empirical interpolation method. SIAM Journal on Scientific Computing, 36(1), A168-A192, 2014.

[17] D. Ryckelynck, F. Vincent, S. Cantournet. Multidimensional a priori hyper-reduction of mechanical models involving internal variables. Computer Methods in Applied Mechanics and Engineering, 225, 28-43, 2012.

[18] D. Amsallem, M.J. Zahr, Y. Choi, C. Farhat. Design optimization using hyper-reduced-order models. Structural and Multidisciplinary Optimization, 51(4), 919-940, 2015.

[19] K. Willcox. Unsteady flow sensing and estimation via the gappy Proper Orthogonal Decomposition. Comput. Fluids, 35(2), 208226, 2006.

[20] P. Benner, S. Gugercin, K. Willcox. A survey of projection-based model reduction methods for parametric dynamical systems. SIAM Review, 57(4), 483-531, 2015.

[21] 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.

[22] A.T. Patera, G. Rozza. Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations. Copyright MIT (2006-2015). MIT Pappalardo Monographs in Mechanical Engineering, online at http://augustine.mit.edu, 2007.

[23] M. Drohmann, B. Haasdonk, M. Ohlberger. Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation. SIAM Journal on Scientific Computing, 34(2), A937-A969, 2012.

[24] F. Fritzen, M. Leuschner. Reduced basis hybrid computational homogenization based on a mixed incremental formulation. Computer Methods in Applied Mechanics and Engineering, 260, 143-154, 2013.

[25] G. Rozza, D.B.P. Huynh, A.T. Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Archives of Computational Methods in Engineering, 15(3), 229-275, 2008.
[26] G. Rozza. Fundamentals of reduced basis method for problems governed by parametrized PDEs and applications. In CISM Lectures notes "Sepa-rated Representation and PGD based model reduction: fundamentals and applications", P. Ladeveze and F. Chinesta Eds., Springer Vienna, 2014.

[27] J. Hesthaven, G. Rozza, B. Stamm. Certified reduced basis methods for parametrized partial differential equations, Springer Briefs in Mathematics, Springer, 2015.

[28] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Journal of Non-Newtonian Fluid Mechanics, 139, 153-176, 2006.

[29] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representation. Journal of Non-Newtonian Fluid Mechanics, 144, 98-121, 2007.

[30] A. Falco, A. Nouy. A Proper Generalized Decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach. J. Math. Anal. Appl., 376, 469-480, 2011.

[31] P. Ladevèze, The large time increment method for the analysis of structures with nonlinear constitutive relation described by internal variables, Comptes Rendus Académie des Sciences Paris, 309, 1095-1099, 1989.

[32] L. Boucinha, A. Gravouil, A. Ammar. Spacetime proper generalized decompositions for the resolution of transient elastodynamic models. Computer Methods in Applied Mechanics and Engineering, 255, 67-88, 2014.

[33] F. Chinesta, A. Ammar, A. Leygue, R. Keunings. An overview of the Proper Generalized Decomposition with applications in computational rheology. Journal of Non Newtonian Fluid Mechanics, 166, 578-592, 2011.

[34] W. Hackbusch. Tensor spaces and numerical tensor calculus. Springer Series in Computational Mathematics. Springer-Verlag, Berlin-Heidelberg, 1st edition, 2012.

[35] P. Lavedèze, L. Chamoin. On the verification of model reduction methods based on the Proper Generalized Decomposition. Comput. Meth. Appl. Mech. Engrg., 200(23-24), 203-22047, 2011.

[36] A. Nouy. A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations. Computer Methods in Applied Mechanics and Engineering, 196, 4521-4537, 2007.
[37] B. Bognet, A. Leygue, F. Chinesta. Separated representations of 3D elastic solutions in shell geometries. Advanced Modelling and Simulation in Engineering Sciences, 1(1), 1-4, 2014.

[38] F. Chinesta, A. Leygue, F. Bordeu, J.V. Aguado, E. Cueto, D. Gonzalez, I. Alfaro, A. Ammar, A. Huerta, PGD based computational vademecum for efficient design, optimization and control. Archives of Computational Methods in Engineering, 20(1), 31-59, 2013.

[39] F. Chinesta, A. Leygue, B. Bognet, Ch. Ghnatios, F. Poulhaon, F. Bordeu, A. Barasinski, A. Poitou, S. Chatel, S. Maison-Le-Poec. First Steps Towards an Advanced Simulation of Composites Manufacturing by Automated Tape Placement. International Journal of Material Forming, 7(1), 81-92, 2014.

[40] D. Borzacchiello, J.V. Aguado, F. Chinesta. Reduced order modelling for efficient process optimisation of a hot-wall chemical vapour deposition reactor. International Journal of Numerical Methods for Heat and Fluid Flow. In press.

[41] M. R. Kamal, S. Sourour, Kinetics and thermal characterization of thermoset cure. Polymer Engineering & Science, 13(1), 59-64, 1973.

[42] B. Gebart. Permeability of unidirectional reinforcements for RTM. Journal of Composite Materials, 26(8), 1100-1133, 1992.

[43] M. Barrault, Y. Maday, N.C. Nguyen, A.T. Patera. An "empirical interpolation" method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathématique, 339(9), 667-672, 2004.

[44] C. Germoso, J.V. Aguado, A. Fraile, E. Alarcón, F. Chinesta. Efficient PGD-based dynamic calculation of non-linear soil behavior. C.R. Mécanique, 344(1), 24-41, 2016.
## About the authors

**Domenico Borzacchiello, Jose V. Aguado and Francisco Chinesta**

Institut de Calcul Intensif (ICI) at Ecole Centrale de Nantes
1 rue de la Noë, BP 92101, 44321 Nantes cedex 3, France
e-mail: {domenico.borzacchiello, jose.aguado-lopez, francisco.chinesta}@ec-nantes.fr
web: http://ici.ec-nantes.fr

**Chady Ghnatios**

Notre Dame University-Louaize
P.O. Box: 72, Zouk Mikael, Zouk Mosbeh, Lebanon
e-mail: cghnatios@ndu.edu.lb

**François Lebel and Ram Upadhyay**

Composite Technologies, GE Global Research
One research circle, Niskayuna, NY 12309
e-mail: {francois,upadhyay}@ge.com

**Christophe Binetruy**

Institut de Recherche en Génie Civil et Mécanique (GeM UMR CNRS 6183),
Ecole Centrale de Nantes. 1 rue de la Noë, BP 92101, F-44321 Nantes cedex 3,
France.
e-mail: christophe.binetruy@ec-nantes.fr