Bending behaviour of thermoplastic composites in melt: A data-driven approach

The resistance to bending of continuous fibre-reinforced thermoplastic composites at processing temperature is an important predictor of wrinkle formation during the stamp forming process. This resistance can be quantified with experiments and approximated with models. However, current models are validated with one, or at most two, thermoplastic composites. Hence, it is not known whether these models generalise to other thermoplastic composites. This data-driven study proposes a rate- and temperature-dependent bending model derived from the experimental data of 10 very different thermoplastic composites. The multivariate data structure was analysed with a Tucker decomposition, together with bootstrapping and cross-validation, to reveal the important gen- eralisable trends, and discard irrelevant features. The most suitable Tucker model could be parametrised to yield a four-parameter model, which approximates the deformation, rate, and temperature dependence, respectively using a linear, a power law, and an Arrhenius-type relation. The four-parameter model was validated with experimental data from four other thermoplastic composite materials. It can describe the measured rate- and temperature-dependent bending behaviour of a wide range of thermoplastic composites


Introduction
With their efforts to increase vehicular energy efficiency, the aerospace and automotive industry are driving the demand for lightweight materials such as continuous fibre reinforced polymers [1]. The challenge is to produce parts from such materials in high volumes. Thermoforming of thermoplastic composites (TPCs) is a promising automated process that is suited for high volume production. However, the forming process can introduce defects, such as wrinkles, waviness, tears and fibre buckling. Wrinkling has a severely negative influence on the part's performance [2], and is one of the most common defects observed [3]. Since the out-of-plane bending behaviour has been shown to be particularly important in predicting the size and shape of wrinkles [4], research efforts have focussed on developing an experimental method to characterise and model the out-of-plane bending behaviour.
The Kawabata Evaluation System for Fabrics-inspired [5,6] bending setup developed at the University of Twente [7], is one of the available methods. It is represented schematically in Fig. 1a. The method quantifies the reaction moment M c of a TPC specimen to an imposed bending deformation, in the form of a rotation α, in a controlled environment.
Typically, a full 3 × 3 factorial design is applied (Fig. 1b), to be able to determine the rate ω and temperature T dependence of the bending behaviour.
To model the bending behaviour measured with the Kawabatainspired bending setup, two approaches can be recognised. First, a model can be based on the micromechanical considerations. Sachs [7] derived a viscoelastic model to describe the bending behaviour of unidirectional (UD) materials, by considering the specimen as a layered structure of elastic fibre filaments and a viscous polymer matrix. The resulting physical model could be formulated as a non-linear diffusion equation in terms of the rotation angle. By adopting a shear thinning power law description for the matrix viscosity, this three-parameter constitutive model can describe the rate dependence of the bending experiments at a constant temperature [8], eg. Fig. 1c.
Second, a rheological model, i.e. a combination of a number of elastic spring and viscous dashpot elements in series or in parallel [9], can be fitted to the experimental results. For example, Ropers et al. [3] found a seven element non-linear viscoelastic model with a total of 18 parameters for each experimental temperature. Their model was able to represent the bending experiments at the three rates, as well as the subsequent relaxation. Dörr et al. [10] compared and evaluated the performance of a two-element Kelvin-Voigt model with five parameters, and a five-element generalised Maxwell model with nine parameters. Both viscoelastic models were able to describe the rate dependence at a constant temperature. As expected, the fit of the generalised Maxwell model proves a better match with the experiments, due to the larger number of fitting parameters.
The previously discussed models can be used to represent the moment-rotation angle curves of the bending experiments performed at different rates at a constant temperature, for instance the measurements presented in Fig. 1c. Furthermore, various approaches were presented to also describe the observed temperature dependence, for instance, by applying the principle of time-temperature superposition [3], interpolating between discrete parameters [11], or making the elastic stiffness tensor and the constitutive relation for equivalent flow rate (not its parameters) a function of the temperature [12]. However, the proposed models and approaches were validated with the measurements of one, or at most two different TPCs, since experimental work is expensive and time consuming. Such tailored modelling imposes a great risk due to the lack of generalisation. Therefore, as argued by Ropers et al. [3], the transferability of the proposed models to other TPCs must be verified.
To help address this point, the experimental data of 10 thermoplastic composite materials 1 listed in Table 1 are analysed in this paper. The characterised materials cover a wide spectrum of TPCs, ranging from TPCs designed for automotive purposes, to TPCs for aerospace applications, including fabrics and unidirectional reinforcement. The aim of the analysis is to derive a model that is able to describe the general relation between the measured moment and the three experimental variables: rotation angle, rotation rate, and temperature, for all these TPC materials. The proposed model does not include geometry dependencies, except for the specimen thickness.
The analysis consists of two parts. First, multivariate data analysis methods are applied to identify the relations between the moment and each of the variables common to all materials. At this stage the most suitable model is selected. This model captures the most important trends, but discards uncommon and irrelevant features. Second, this model is parametrised, and further simplified. Then, the final dataderived model is fitted to the original data, and validated with new data. Finally, a physical interpretation of the data-derived model is discussed.

Experimental data
On the Kawabata-inspired bending setup (Fig. 1a), a specimen is loosely placed between two clamps allowing for a sliding movement during the experiment. The setup is heated in the nitrogen rich thermal chamber (CTD 450) of the Anton Paar rheometer (MCR 501) to well over the melting temperature of the TPC matrix material. After 3 minthree minutes, a rotation α is applied, from 0 • up to 70 • , while the moment M c is measured. Typical measurement curves are plotted in Fig. 1c. The rotation rate ω is constant from 10 • to 65 • . More details on the experiment can be found in Ref. [8].
The bending behaviour of a TPC m is characterised by measuring the moment for three rates, 0.6 • /s, 6 • /s and 60 • /s, at three temperatures T m − 20 • C, T m , and T m + 20 • C, where T m is the recommended processing temperature ( Table 1). The grid of nine experiments is shown in Fig. 1b. Normally, each experiment is repeated on three different specimens, in practice however, between two and five valid curves are available. For this research, as mentioned earlier, the characterisation data of 10 different TPC materials are analysed.

Data pre-processing
In the pre-processing stage, the experimental data are scaled, cropped, and resampled. The moment is divided by the number of plies of the specimen (Table 1), to obtain the measured moment for one ply [8]. The measurement curves are then cropped to remove the acceleration and deceleration phases, hence only the measurements for rotation angles from 10 • to 65 • are included. Moreover, the data are resampled equidistantly in intervals of 1 • by means of linear interpolation. Finally, any missing values are masked in such a way that these can later be automatically estimated by the algorithms used.
The data array , in which the processed data are stored, is constructed to include the material type, deformation (rotation angle), rotation rate, temperature, and specimen number, each as an independent dimension. Therefore, any data point x mdrtp in the dataset is uniquely specified by its indices: m = 1, …, M, d = 1, …, D, r = 1, …, R, t = 1,…,T, and p = 1,…,P mrt , where M = 10, D = 56, R = 3, T = 3, and P mrt = 2, …, 5. In this study, the vector with the 56 interpolated measurements of one specimen will be written as x mrtp .

Exploratory data analysis
The Tucker decomposition [13] is applied in conjunction with bootstrapping [14] and cross-validation [15] to model the relationship between the measured moment and the four controlled sources of variation: material, deformation, rotation rate, and temperature, while taking the random variation between repeated experiments into account.

Tucker decomposition
The Tucker decomposition is a generalisation of the principal component analysis to n-dimensional arrays [16], and is used to decompose an n-dimensional data array 2 into n factor matrices A and an n-dimensional core array , based on the mutual linear correlations in the dataset. For this analysis, four-dimensional data arrays are Fig. 1. a) Schematic drawing of the bending setup. Reprinted from Sachs and Akkerman [8], with permission from Elsevier. A clamp rotation α is applied, while the moment M c is measured. The distance between the pivot and the clamp l is 7.5 mm. b) Experiments are performed at nine temperature-rate combinations. T m is the recommended processing temperature of the TPC. c) Typical measurement curves at a constant temperature, and three constant rotation rates. 1 The warp and weft direction of the 4-1 biased plain weave glass fibre PA66 are counted as separate materials. 2 In the statistical literature multidimensional data arrays are often referred to as tensors, but these are not constrained as tensors in continuum mechanics are.
considered ∈ R M×D×R×T . Their four dimensions correspond to the first four dimensions of the dataset, hence, the size of each data array is 10 × 56 × 3 × 3. For such an array the Tucker decomposition reads: in which ̃i s the estimate of the input array, × N is the n-mode 3 tensor product [17,18]. In index notation, the decomposition in Equation (1) is written as: The columns of the factor matrices contain the principal components for each dimension in order of importance. The first few principal components capture the significant trends in the data, whereas the last components contain the noise. Hence, a truncated Tucker model, where J ≤ M, C ≤ D, Q ≤ R, and S ≤ T, can often give a good approximation of the data array. Such a truncated model is denoted as a Tucker(J, C, Q, S) model.
The Tensorly implementation [19] used in this research decomposes the data array with the help of Higher Order Orthogonal Iteration, as described by Kolda and Bader [16]. The algorithm can estimate missing values and is initialised with a random core array and random factor matrices.

Bootstrapping
The bootstrap method [14] was used in this analysis as a means to estimate the prediction error of the Tucker models, as well as the uncertainty of the principal components. A total of 1000 bootstrap samples 4 {b} , with b = 1, …, 1000, are generated from the original dataset . To create a bootstrap sample, each of the 90 vectors y {b} mrt of {b} is filled by selecting one of the P mrt available repeated specimen x mrtp at random from the matrix X mrt with replacement. Hence, the vector with the measurements of a specimen remains intact. This retains the relation between subsequent measurements along the deformation dimension, while the relation between specimens along the material, rate and temperature dimensions will vary. The data array complementary to {b} will be denoted as {b} .
Every bootstrap sample is decomposed by use of Equation (1) to its core array and factor matrices. For example, the components obtained by decomposing the 1000 bootstrap samples of the Tucker (5,5,3,3) model are plotted in Fig. 2. Each row of graphs presents a factor matrix, with each of its components (i.e. columns) in one graph. Because the bootstrap samples differ slightly, variation in the components also exists. The first components have little variation, and are almost deterministic. These are the significant trends common to all observations. In contrast, the variation is clearly visible for the latter, less important and noisier, components.

Model selection
Cross-validation is applied to determine the appropriate truncated Tucker(J, C, Q, S) model, based on the out-of-sample predictive performance. In this context, every truncation is a separate model, and the sum of the components (J + C + Q + S) will be referred to as the model size. The predictive performance is determined by dividing the dataset into two sets, a train set and a test set. In this study, the selected specimen in a bootstrap sample {b} are the training data, while the remaining repeated experiments can be used as test data {b} . The dataset for training is decomposed to get the estimate ̃{ b} . The error of the insample prediction, or training error, is the mean squared error (MSE) between the estimate ̃{ b} and the training data {b} . The out-of-sample prediction, or test error, is the MSE between ̃{ b} and the test data {b} . Training and testing are repeated on all bootstrap samples. Therefore, a distribution of the train and test error is obtained per model. The sum of components is set to range from 6 to 20, while the number of components per dimension is limited to J ≤ 10, C ≤ 15, Q ≤ 3, and S ≤ 3.
The train and test error distributions of the best performing model for each model size are plotted in Fig. 3. The train and test error quickly drop as the total number of components increases. The test MSE levels off at around a model size of 10. The Tucker(5, 3, 2, 2) model, with a model size of 12, emerges as the best performing model. However, there is no significant difference with respect to the smaller Tucker(4, 2, 2, 2) model. Moreover, the models around the best performing model all suggest that the general relation between the moment and the deformation, rate, and temperature can be captured with two components. These components have an effect on the level, and the slope of the moment relation [20]. Furthermore, only four, or at most five, material Table 1 The overview of the reinforcement type, fibre material, and matrix polymer material of the thermoplastic composites tested, as well as the index number m and reference used in this paper, the number of plies per specimen, ply thickness h, and recommended processing temperature T m . Furthermore, the specimen width W is 25 mm and its length L is 35 mm. The matrix polymer materials are abbreviated as follows:   3 The n th dimension of an array is often referred to as the n th mode. 4 The number of bootstrap samples was chosen based on the convergence of the distributions of the coefficients of the factor matrices, and the significant figures reported.
coefficients are necessary to describe the bending behaviour of a TPC material. Therefore, the Tucker(4, 2, 2, 2) model was selected as the most suitable model. Its components are coloured blue in Fig. 2, and are indistinguishable from those of the Tucker(5, 5, 3, 3) model.

Parametrisation
The deformation-, rate-, and temperature-related components of the chosen Tucker(4, 2, 2, 2) model can be parametrised to facilitate the model's interpretation and application. The first two deformation components could be approximated with linear equations. A power law and an Arrhenius-type relation were applied to describe the non-linear rate and temperature dependence respectively [8,9]. The parametrised components, denoted with the hat ^, can be written as follows: â r rq = β r 0q + β r 1q ω n r r = 1, 2, 3, q = 1, 2, â t ts = β t 0s + β t 1s e Eη RTa t t = 1, 2, 3, s = 1, 2, with the regression coefficients β, the flow behaviour index n, the activation energy E η , and the gas constant R. At this stage, the average processing temperature T a is used, where T a 2 = 593.15 K. Note that the flow behaviour index and the activation energy are shared by both components, and that the coefficients are material independent. The fit, i.e. the dashed light blue line in Fig. 2, shows solid alignment with the specified components. The rate and temperature estimates show more uncertainty, i.e. the blue shaded area, as there are only three points to fit.
Since there were four material coefficients, a M m1 to a M m4 , in Equation (6), only four of the eight material-dependent coefficients can be considered independent. Therefore, the number of material-dependent parameters in Equation (7) can be reduced. In the following steps Equation (7) is reduced to a equation with four material-dependent parameters based on the following two observations. First, for most materials there is a strong correlation (mean 0.97, min. 0.75, max. 1.00) between the intercept and the slope. This relation between the intercept and the slope could be approximated with a proportional relationship. Hence, Equation (7) can be approximated with: in which the parameter ρ m scales the slope relative to the intercept.
Second, the rate and temperature dependence of the intercept is captured by the four terms, with one flow behaviour index n, one activation energy E η , and an average processing temperature T at for all materials. This dependence can be approximated by a power law--Arrhenius-type equation per material with its actual processing temperature T m and a material specific flow behaviour index n m and activation energy E ηm . With this knowledge, Equation (8) can be written as: where the proportionality factor k m , flow behaviour index n m , activation energy E ηm , and slope factor ρ m constitute the four parameters per material m. Again, it is assumed that the moment scales linearly with the specimen thickness H m [8]. Finally, following Equation (9)

Parameter inference
Bayesian inference [21] is applied to determine the distributions of the four material parameters in Equation (10) for each of the 10 TPC given the experimental data . Further, the parameters of four additional TPC materials-W5HS-CF-PPS, W22TW-GF-PA6 R and UD-CF-PA6 R from Ropers et al. [3,11], and UD-CF-PA6 S used by Sachs and Akkerman [8]-have been estimated as well. The implemented model aims to estimate the mean of each repeated experiment. Although this model does not model the data generating process of the samples, it does indicate where the data for the samples may be expected. The model is defined as follows: [model] (12) σ ∼ Normal(0, 1) [prior] n ∼ Beta (2,2) [prior] where the estimate of the moment M est (Equation (11)) is assumed to be normally distributed with a mean μ (Equation (12)) and a standard deviation e σ . Equation (12) is a reparametrisation of Equation (10) that improves sampling. To obtain Equation (12), the deformation, rate, and temperature relations are normalised so that their relative contribution is one at α 28 = 37π/180, ω 2 = 6π/180, and T 2 = T m respectively, i.e.
the midpoint of the experimental design. Additionally, the proportionality factor k and slope factor ρ have been log-transformed and scaled to compensate for the normalisation. These transformed and scaled parameters are denoted with k and ρ. This reparametrisation has the added advantage that it facilitates the comparison between the proportionality factors of different TPCs. Because, the transformed proportionality factor k is independent of the rate and temperature range considered. In contrast, the proportionality factor k in Equation (10), depends on the processing temperature, activation energy, and flow behaviour index of the TPC, since The standard deviation (Equation (13)) and four parameters (Equations (14) to (17)) are given informative priors. The inference was performed with the No U-Turn Hamiltonian Monte Carlo [22] sampler of the PyMC3 package [23]. For details on the implementation of the model and sampler settings, please refer to the supplementary materials.

Results
The estimate of the four-parameter model captures the general trends, with respect to the deformation, rate, and temperature, observed for all the TPCs materials. For instance, the estimate for the measurements of the UD-CF-PA6 and UD-GF-PA6 composites are plotted in Fig. 4, top and bottom row, respectively. Figures with the estimates for the other TPCs, as well as the means and covariance matrix of the parameters, can be found in the supplementary materials.
The plots show the data measured at the nine experimental settings in grey. The blue lines represent the estimates for the mean of the experiments, and the light blue shaded area is the 95% compatibility interval. The dark blue shaded area represents the compatibility interval for the estimate mean μ, i.e. 95 % of the curves generated by the model, with random samples drawn from the parameter distributions presented in the supplementary materials, will lie in this shaded area. Note that the four-parameter model does not fully capture uncommon features, such as, pronounced curvatures (Fig. 4a), or slopes independent from the rate or temperature (Fig. 4d-f).
The distributions for the four material parameters are presented in Fig. 5. A clear distinction can be made between the parameters of each material. The proportionality factor k scales the estimate. TPCs on the left have, on average, a lower resistance to deformation than those on the right. The remaining parameters ρ, n and E η , define the sensitivity of a TPC to deformation, rate, and temperature changes respectively. TPCs on the left are less sensitive to change, compared to TPCs on the right. Although there is significant variation between the materials for these parameters, their effect on the moment is relatively small compared to the effect of the proportionality factor. For instance, the effect of the proportionality factor on the estimated moment, ek, was found to range from e − 3.05 = 0.047 4 (UD-GF-PA6) to e 1.65 = 5.21 (UD-CF-PEEK), which is a difference of a factor 110. In contrast, the effect the lowest and highest parameter values of ρ, n and E η have on the moment differ by a factor of at most 16, 2, and 2, respectively.

Discussion
An analysis of the experimental bending data was performed to derive a descriptive model that generalises to a wide variety of materials, and describes the relation between the measured moment and the rotation angle, rotation rate, and temperature. The appropriate model was determined in the first part of the analysis by applying the Tucker decomposition on the gathered experimental data in combination with bootstrapping and cross-validation. The Tucker(4, 2, 2, 2) model was found to be the most suitable model. In the second part, the chosen Tucker model was parametrised and simplified, to yield an equation with four parameters. A linear relation was found between the moment and the rotation angle. Interestingly, the slope of this linear relation is rate-and temperature-dependent and is proportional to the intercept. The rate dependence could be described by a power law type relation, and an Arrhenius-type equation was applied to describe the temperature dependence. The model was shown to describe the bending behaviour of the 10 thermoplastic composites from which the model was derived, as well as four additional TPCs. Since, the presented model is solely based on experimental data, it can only model the bending behaviour of this specific setup and specimen geometry for rotation angles between 10 • and 65 • , at constant rotation rate and temperature. Note that the exploratory analysis method described may be extended to describe the full measurement curve from 0 • to 70 • . For instance, by applying a continuous piecewise differentiable function to describe the deformation components.
The four parameters-k, n, E η , and ρ-of the final model, as presented in Equation (12), each have a distinct influence on the estimate of the moment. The latter three parameters each influence one aspect of the estimate: the flow behaviour index n, and activation energy E η have an effect on the rate, and temperature dependence of the TPC respectively, and the slope factor ρ determines the deformation dependence. On the other hand, the proportionality factor k scales all aspects of the estimate simultaneously. This parameter also captures most of the variation between the different materials, since the normalised deformation, rate, and temperature relations vary little compared to the variation of ek. Therefore, the primary differentiating characteristic between the bending behaviour of TPCs is the mean resistance to deformation, i.e. level of the measurement curves.
The four-parameter model, in the form of Equation (10), resembles the long-term solution of the viscoelastic bending model presented by Sachs and Akkerman [8]: where W, H and L are the width, height and length of the specimen, l is the distance between the pivot point and the clamp, and c, k c , and n are the elastic bending stiffness, flow consistency index, and the flow behaviour index respectively. Both equations (10) and (19) describe a linear relation between the moment and deformation, and a shear thinning power law relation for the rate dependence. Moreover, an Arrhenius-type relation for the flow consistency index k c could be adopted to also make the long-term solution (19) temperaturedependent. In contrast to the earlier solution, here, the slope was found to be significantly rate-and temperature-dependent for most materials, which is consistent with previous observations by other authors [8,11,12]. This may be taken into consideration in future research. The effect of the thickness was assumed to be linear in this study, as derived by Sachs and Akkerman [8] for unidirectional TPCs. The data available to verify this assumption are limited to a few additional measurements for UD-GF-PA6, W22TW-GF-PA66, W41BP-GF-PA66, and W5HS-CF-PPS. The expected moment for these additional measurements at a different thickness could be predicted with Equation (10), using the new thickness and the parameter estimates found for the original thickness. The predictions are presented in the supplementary materials. This suggests that the moment scales linearly with the thickness for both unidirectional and woven TPCs. However, more research is needed to validate this linear relation and to determine the underlying mechanisms.
The data-derived four-parameter model captures the general bending behaviour for the available experimental data. More experiments, or thorough domain knowledge, are needed to justify the choice for more complex models that can capture uncommon features. Furthermore, the required accuracy of the estimate ultimately depends on the sensitivity of the forming simulations, as the research by Dörr et al. [10] illustrates. For the tested geometries in those studies, the five-parameter Kelvin-Voight model, and the nine-parameter generalised Maxwell model produced comparable simulation results. Further research is needed to analyse the sensitivity of the simulation results. Such an analysis would be facilitated by applying the method described in this study to the other deformation mechanisms as well.

Conclusion
A descriptive model was derived that describes the integral response of a thermoplastic composite to bending in forming conditions for the Kawabata-inspired bending setup. This four-parameter model is applicable for a wide range of thermoplastic composite materials. Interestingly, the slope of the linear relation between the moment and rotation angle, was found to be rate-and temperature-dependent and proportional to the intercept. A power law and Arrhenius-type relation were appropriate to capture, respectively, the rate and temperature dependence of the thermoplastic composites. The four parameters of the model can be interpreted as a scale, rate, temperature, and slope parameter. In general, all thermoplastic composites evaluated here show a similar, yet distinct, bending behaviour that can be captured by a simple model. The model can be used to fit new experimental results, to place these results in the broader context of already tested materials, and to predict expected experimental results at rates and temperatures other than those tested. The application of the analysis method to the other deformation mechanisms will be a topic of future research.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.