**1. Introduction**

Sheet Moulding Compounds (SMC) compression moulding is the largest market segmen<sup>t</sup> in the processing of Glass Fibre-Reinforced Plastics (GFRP) in terms of production volume [1]. The benefits of compression moulding include the economical production of near-net-shape components, minimising the need for subsequent assembly steps. The deformation of the SMC preform during compression moulding causes the formation of inhomogeneous and transient flow fields, which in turn cause a change in the orientation of the contained fibres [2–8]. Mechanical properties of the resulting SMC part are dominated by these fibre orientations, which therefore play an important role in the design of SMC parts and setup of the compression moulding process such as definition preform position and size [4,9–11].

Martulli et al. reported a property di fference between specimens cut from carbon fibre-reinforced (53 wt %) vinylester-based SMC with a preferential orientation of 0◦ and 90◦ of 150% for tensile sti ffness, 260% for tensile strength, 120% for compressive sti ffness, and 32% for compressive strength, respectively [11]. For a polyester-based SMC with 30 wt % glass fibre content, Oldenbo et al. determined a di fference in tensile sti ffness for preferential orientation of 0◦ and 90◦ upwards of 25% [10].

Fibre orientation can be predicted by the application of process simulation procedures, on which extensive work has been conducted since the early 1980s [12–16]. Initially, these were based on 2D and 2.5D modelling approaches (the latter considering through-thickness variations in flow and material properties). Lee, Folgar, and Tucker applied the generalised Hele–Shaw model for calculation of the filling of thin-walled structures during the SMC compression moulding process; however, they did initially not consider the influence of temperature on the SMC viscosity [15,17,18]. Barone and Caulk performed experimental analysis on SMC flow and reported boundary e ffects such as slippage of the SMC on the mould surface due to temperature influences, for which a model was subsequently proposed [19,20].

In recent years, these methods were expanded by the development of simulation procedures capable of 3D calculation, of which certain functionalities have been implemented in programs such as Moldex3D and Moldflow [21,22]. Hohberg employed the Coupled Eulerian Lagrangian framework within Abaqus for the calculation of SMC flow and flow-induced deformation of local reinforcements [23]. A research group led by Osswald developed a direct fibre simulation procedure, with which fibre bending and fibre–matrix separation can be calculated in long fibre-reinforced polymers [21]. This was expanded on by Meyer et al. with a direct bundle simulation approach [24].

However, process simulation has remained computationally demanding andtime-consuming [4,9,25]. This limits the potential in part or process optimisation, since finding the "optimal" solution usually necessitates simulating a high number of variable variations, which is especially prevalent when varying only one at a time [25–30]. In the field of SMC compression moulding, advanced optimisation procedures have been presented, which are based on alternative approaches making use of approximations of the complex interactions of influencing parameters (also called metamodels or surrogate models) [25,30,31]. Metamodels provide a "Model of the Model", which may be used to replace computationally expensive simulation models in a wide number of engineering disciplines and have also been used in SMC modelling and process optimisation [25,30,32–34]. Huang et al. used a mesoscale metamodelling approach to accurately predict the sti ffness matrices of chopped carbon fibre SMC from the fibre orientation tensor using individual Kriging models for each element of the sti ffness matrix [34]. Sabiston et al. presented a neural networks (NN)-based procedure for the prediction of preform position and geometry-dependent fibre orientations in a SMC seat back component. With this approach, near instantaneous prediction of fibre orientation is achieved, the caveat being that a high number of data points (3000) is required for training the NN [8].

For determination of the optimal SMC peform placement for a hood scoop part, Ankenman, Bisgaard, and Osswald proposed an iterative optimisation approach based on evaluating simulation results with the response surface methodology [25]. The optimisation had two goals, being to minimise the "fill time tolerance" (standard deviation of the time necessary for filling all nodes of the discretised part geometry, thus describing the filling uniformity) and the uniformness of fibre orientation. However, this procedure necessitated significant human interaction in iteratively defining and evaluating di fferent experimental setups. This was later expanded on by Twu and Lee, who automated the procedure by applying linear regression analysis. However, they noted that this methodology (as well as statistical regression methods or various mathematical approximation theories in general) had major drawbacks, as the number of necessary simulations would increase exponentially with the number of design variables [35]. Alternatively, they proposed employing neural networks (NN), which they subsequently applied in increasing the curing homogeneity of an SMC part with varying thickness by optimising the heating channel location inside the mould [30]. Initially, a start-up search is employed (necessitating 19 curing simulations). The parameter variations and simulation results are used for initial training of a feed-forward NN (FF-NN). Subsequently, the FF-NN is iteratively evaluated and retrained by supplementary curing simulation, finding the optimal design in less than 60 simulations (in comparison with the statistical approach necessitating 729 simulations when using a 3-level quadratic model without domain search) [30].

Alternatively, Kim, Lee, Han, and Vautrin used a genetic algorithm for optimising the preform size and placement. The goal of the optimisation was to minimise the maximum deflection of a symmetric car hood and an arbitrary non-symmetric geometry that resembles a fender [31].

Most of the optimisation procedures mentioned apply initial evaluations of the significance of design factors, excluding non-significant factors, to make the design cases more manageable [25,30]. Additionally, problem-dependent constraint handling techniques are used to rule out non-feasible solutions (e.g., limiting the ranges of the design factors to physical limits such as the mould dimensions). For curing homogeneity optimisation, Twu and Lee applied metamodelling to four design factors [30]. Ankenman, Bisgaard, and Osswald limited the number of design factors to three (charge size, length-to-width ratio, and position relative to one mould edge), while Kim, Lee, Han, and Vautrin employed the penalty function method and a repair algorithm in the handling of four design variables (size and position of the preform in x and y direction) [25,31].

Thus, promising approaches have been shown addressing the presented challenges, in which a wide range of metamodelling procedures are successfully used. However, to the best of the authors' knowledge, no metamodelling procedures that expand on the geometrical freedom of the preform (e.g., non-rectangular preforms) have been presented in the field of SMC processing.

Therefore, the focus of the presented study is the implementation of a metamodelling approach without inherent geometry restrictions, for which an ensemble metamodel comprising multiple FF-NN is proposed. The procedure directly makes use of the spatial geometry discretisation used in 2D and 2.5D process simulation for the description of the part and preform geometry and approximation of the correlation of the preform and resulting mechanical properties. The metamodel is subsequently used in a procedure for the optimisation of the SMC preform. As the design goal, minimisation of the maximum absolute deflection of a plate geometry under bending load is pursued.

Sampling for FF-NN training is conducted by evaluating rectangular preform geometries (which in sum span the totality of the part surface) and the resulting maximum absolute deflection using a coupled process and structural simulation. In the following subsection, methods used for the metamodelling of preform and part behaviour and subsequent metamodel-based preform optimisation are presented.

The presented procedures were implemented in MATLAB R2020a, Mathworks, Natick, MA, USA, when not stated otherwise.

#### **2. Materials and Methods**

The procedures were applied to a plate geometry exhibiting a cantilever load case, which is shown in Figure 1a. The goal of the optimisation of the preform was the minimisation of plate deflection. As this paper concentrates on procedure development, this fairly simple geometry is chosen.

**Figure 1.** Minimisation of the plate's maximum absolute deflection: Plate geometry, boundary conditions, and applied load (**a**). Discretisation of plate geometry with 1200 S3 elements (**b**).

Discretisation of the geometry into a structured shell mesh was conducted in Abaqus 2020, Dassault Systèmes, Vélizy-Villacoublay, France, using 1200 S3 elements (Figure 1b). This shell element type was chosen based on compatibility with the process simulation procedure, and subsequently, it was also used in structural simulation and metamodelling.

### *2.1. Explicit Modelling*

Calculation of the fibre orientation probability distribution function (FOD) ψ(*p*) resulting from SMC flow during compression moulding was conducted using Express 6.0, M-Base Engineering + Software GmbH, Aachen, Germany. This software was developed in close collaboration with the Institute for Plastics Processing (IKV) (Aachen, Germany) and it has been used for a 2.5D process simulation of thermoplastic and thermoset compression moulding alike [36–38]. It is based on the control volume approach as described by Osswald, which has been shown to accurately predict the filling pattern in compression moulding of thin geometries under the assumption of planar flow [12]. The governing equations implemented in the software and used in this paper are described briefly. The material data used are shown in Section 2.2.

To predict the static pressure distribution p within the mould cavity during the compression moulding of SMC, Folgar and Tuckers' application of the generalised Hele–Shaw model is used [12,15,17]:

$$\frac{\partial}{\partial \mathbf{x}} \left( S \frac{\partial p}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{y}} \left( S \frac{\partial p}{\partial y} \right) - \dot{h} = 0 \tag{1}$$

The flow conductivity *S* is derived from the flow gap height *h* and the shear and temperature-dependent viscosity η:

$$S = \frac{h^3}{12\eta} \tag{2}$$

.

From the pressure distribution, gap-wise average velocities *U* and *V* are derived, from which the fibre orientation probability distribution function (FOD) is subsequently calculated [12,15]:

$$\overline{M} = -\frac{\mathcal{S}}{h} \frac{\partial p}{\partial \mathbf{x}} \tag{3}$$

$$\overline{V} = -\frac{S}{h} \frac{\partial p}{\partial y} \tag{4}$$

Anisotropy of temperature (and therefore viscosity) is taken into account by simultaneously simulating the filling of the geometry with five shell geometries (assumed as each having one-fifth of the thickness of the plate geometry and being stacked on top of each other), for which the momentary temperature change due to conductive heat transfer is considered individually.

From the layer-wise average velocities, fibre orientations are determined for each layer. Je ffery developed a procedure that enables the prediction of change in orientation of a single ellipsoidal particle due to the flow of a surrounding fluid [39]. For the calculation of the FOD of fibre-reinforced materials, Folgar and Tucker supplemented this model with a phenomenological di ffusion term, with which fibre–fibre interactions are taken into account by the fibre interaction coe fficient *CI* [13]:

$$\frac{\partial \psi}{\partial t} = \text{Cl} \frac{\partial^2 \psi}{\partial \Phi^2} - \frac{\partial}{\partial \Phi} (\psi (-\sin \Phi \cos \Phi \frac{\partial v\_x}{\partial x} - \sin^2 \Phi \frac{\partial v\_x}{\partial y} + \cos^2 \Phi \frac{\partial v\_x}{\partial x} + \sin \Phi \cos \Phi \frac{\partial v\_y}{\partial y})) \tag{5}$$

In the past, extensive work has been conducted in enhancing this orientation model (e.g., by including fibre–matrix interaction), which has resulted in the development of new models such as Anisotropic Rotary Di ffusion (ARD), Reduced Strain Closure (RSC), ARD-RSC, and an improved ARD model combined with the Retarding Principal Rate model (iARD-RPR) [40–45].

In this work, the Folgar–Tucker model is used as recent investigations conducted by Li, Chen et al. continue to show good agreemen<sup>t</sup> with the experimentally determined material properties of SMC and the model having been used in prior work on SMC part thickness optimisation by Kim et al. [22,46].

FOD were transferred to MATLAB for calculation of the mechanical properties of each element in each layer. The mechanical properties are calculated by methods described by Advani and Tucker, which are suitable for use in thin-walled compression moulding (thus assuming planar fibre orientation) and subsequently implemented and validated by Oldembo et al. for SMC materials [5,10]. The calculation of mechanical properties is divided into three successive steps. Initially, the fourth-order stiffness tensor *Cii* (Voigt notation) is calculated under an assumption of unidirectional fibre orientation in the SMC using the governing equations of Halpin and Tsai [47]. Fibre orientation tensors of second order *aij* and fourth order *aijkl* are subsequently derived from the FODs [5,48]:

$$a\_{ij} = \oint p\_i p\_j \psi(p) dp \tag{6}$$

$$a\_{i\bar{j}kl} = \oint p\_i p\_j p\_k p\_l \psi(p) dp. \tag{7}$$

Then, the FOD-dependent stiffness tensors T of each element are calculated by orientation tensor averaging [5]. Scalars *Bi* are derived from the unidirectional stiffness tensor (governing equations in Appendix A), where δ is the Kronecker-delta [5]:

$$T\_{ijkl} = B\_1(a\_{jkl}) + B\_2(a\_{i\uparrow}\delta\_{kl} + a\_{k\ell}\delta\_{l\uparrow}) + B\jmath\{a\_{k\uparrow}\delta\_{jl} + a\_{l\ell}\delta\_{jk} + a\_{j\ell}\delta\_{kl} + a\_{j\ell}\delta\_{l\uparrow}\} + B\jmath\{\delta\_{l\uparrow}\delta\_{kl}\} + B\jmath\{\delta\_{k\uparrow}\delta\_{jl} + \delta\_{l\downarrow}\delta\_{jk}\} \tag{8}$$

For final simulation of part deformation behaviour, the resulting stiffness tensor components were exported to an Abaqus 2020 input-file (.inp) as individually defined materials for each element (thus creating 6000 individual materials) using an automated script. The 5 layers were treated as plies of a composite shell section, which was also created using the script. As fibre orientations are provided in local coordinate systems in the control volume approach, these were also supplied to the .inp file for each element [12].
