Surrogate Model Based Analysis of Inter-Ply Shear Stress in Fiber Reinforced Thermoplastic Composite Sheet Press Forming

: The anisotropic nature of ﬁber reinforced composite materials causes great challenges in predicting the inter-ply shear stress during forming. The complexity of understanding the functional dependency of inter-ply shear stress on multiple forming parameters such as blank temperature, pressure load, inter-ply slippage, and the relative ﬁber orientation angle of adjacent plies further limits the effort to produce a defect-free composite structure. Performing real experiments for various combinations of the mentioned parameters is both time consuming and economically costly. To overcome these difﬁculties, a surrogate-based analysis of inter-ply shear stress is proposed in this study. Based on the ranges of the forming parameters, computer experiments were performed. Using these experimental data, a radial basis function (RBF) based surrogate model that mimics inter-ply shear stress during composite press forming was developed. The ﬁdelity of this model was checked with test data and found to be over 98% efﬁcient.


Introduction
As a result of their superiority in stiffness, impact and corrosion resistivities, as well as being lightweight, reinforced thermoplastic composite materials gain profound advantages over metallic alloys especially in aerospace industries. Unlike their metallic counterparts, reinforced composite materials can be tailored to meet a desired structural performance taking fiber orientation or ply stacking sequences as design variables. Marco Montemurro et al. [1] discussed this possibility with the least-weight design objective constrained to mechanical, geometrical, and technological requirements.
Nevertheless, these privileges come with defect-free forming challenges. Regardless of the advancements in forming technology, manufacturing defect-free composite structures is not yet realized, and hence, the advantages could not be fully exercised. James M. Finley et al. [2] discussed the influence of defects on mechanical performance. Based on their virtual testing framework, defects such as fiber misalignment significantly influence the mechanical response of the composite structure. Various physical experiments and computer simulations have been done over the years [3][4][5][6], though the objective of suppressing forming defects still remains unachieved. Panding et al. [7] studied the influence of micro-voids on the mechanical properties of plain weave carbon/epoxy composite. They reported that the stiffness of the defective composite laminate reduces compared to its theoretical value.
Tim Frerich et al. [8] studied the limitations and challenges of manufacturing impact resistant composite structures. To enhance the impact resistance of composite structures, interleaf layers are inserted between adjacent composite plies. This induces forming defects in composite structures.
A comprehensive review on the manufacturing defects and their effects on the mechanical performance of fiber reinforced composite materials was given by Mahoor Mehdikhani et al. [9]. In their report, it was indicated that, despite many years of research on the formation and effects of defects in composite structures, no mature predictive tools are realized yet.
As long as manufacturing a completely defect-free composite structure is unlikely, to prevent a catastrophic in-service failure of this structure, a post-manufacturing defect detection mechanism has to be developed. In this realm, detection mechanisms have been developed over the years. A case study on wrinkle detection mechanisms reported by Beatriz Larranaga-Valsero et al. [10] depicted both the challenges and essentiality of commonly used detection mechanisms. A feasibility study on a damage detecting fabric sensor for draped composite structures was presented by Jin-Hun Bae et al. [11]. Jun et al. [12] developed an analytical model and analyzed the influence of out-of-plane wrinkle defects on the mechanical properties of fiber reinforced composite laminates. Their report showed that out-of-plane wrinkle defects drastically affect the shear moduli and Poisson's ratios of composite laminates.
Fiber reinforced composite structures can be manufactured through different methods that range from tedious hand-layup to fast and automated fiber placing (AFP) technology. Even with advanced manufacturing technology like AFP, forming defects exist [13,14]. Another commonly known forming technique is press forming, which is used for fast composite structure production especially for aircraft components. Sikken stiffener type, landing gear, door, and rear pressure bulkhead aircraft components are some examples of components manufactured using press forming. Mitsuhiro Okayasu et al. [15], using their hot-stamping system, reported the influence of forming temperature on the incursion of defects. During press forming, adjacent plies of composite laminate have to slip over each other to conform to a mold geometry. Resistance to this slippage results in fiber misalignment, buckling, and wrinkling, thereby degrading the advantages of composite materials. Through experimental and numerical investigations, De Lucas [16] concluded that defects like fiber buckling and wrinkling are caused by inter-ply shear stress induced during shape forming.
To minimize this inter-ply shear stress, adjacent plies have to easily slip over one another. Research findings depict that during composite forming, inter-ply slippage relaxes inter-ply shear stress, thereby reducing fiber buckling and wrinkling [17,18]. Vanclooster et al. [19] conducted experimental analysis about the effects of blank temperature, pressure load, and ply pullout speed on the inter-ply shear properties of a glass/PPcomposite blank and drew the conclusion that shear stress increases for lower temperatures and reduces for decreasing ply pullout speed. Sjolander et al. [20] investigated plausible causes of wrinkling defects during forming. According to their report, the relative fiber orientation angle of adjacent plies has a remarkable effect on wrinkling development. Lightfoot et al. [21] reported the mechanism through which wrinkle is incurred during composite forming. In their report, it was shown that the incursion of wrinkles is also related to in-plane fiber misalignment. Erland et al. [22] presented a technique for determining the forming parameters such as temperature, forming speed, and pressure required to minimize inter-ply shear resistance.
Haanappel et al. [23] conducted a comparative study on uni-directional polyetheretherketone (UD/PEEK) and 8-harness satin polyphenylene sulfide (8HS/PPS) composite materials and concluded that the materials with higher inter-ply shear stress cause severe wrinkling during doubly-curved structure forming. However, due to multiple forming parameters such as the relative fiber orientation angle of adjacent plies, blank temperature, pressure load, and inter-ply slip velocity, it is very complex to determine how and to what extent inter-ply shear stress should be reduced to achieve defect-free products. Scherer et al. [24], in their one-dimensional ply pullout experiment, listed out forming speed, blank temperature, and pressure load as the parameters influencing inter-ply slippage. Hallander et al. [25], in an attempt to achieve wrinkle-free forming of a composite laminate, reported that the relative fiber orientation angle of adjacent plies has a significant influence on the formability of composite structure.
Performing physical experiment for the ranges of these forming parameters is economically costly and time consuming. To ease this costly task, developing a surrogate model that emulates the functional dependency of inter-ply shear stress on forming parameters is vital. Implementing surrogate model based analysis in composite materials becomes a disparate alternative to tedious physical experiments. This paradigm shift is apparent in many of recently published papers [26][27][28][29]. Tifkitsis et al. [26] applied a surrogate model to the cure process of composite laminate so as to remarkably reduce the computational time. Ali et al. [27] also implemented a surrogate model to approximate the mechanical response of composite laminate, thereby demonstrating the computational efficiency and accuracy of the surrogate model. Bassam et al. [28] employed a surrogate model to the meso-scale mechanical response of a composite structure containing wrinkle defects. In an attempt to determine the performance of defective composite structures and their failure modes, Davidson et al. [29] utilized a surrogate model.
Hence, in this study, the mathematical model of inter-ply shear stress, as a function of the aforementioned forming parameters, was formulated and analyzed. The formulation of this mathematical model is based on the surrogate modeling technique. The surrogate model construction and validation procedure taken in this study are as follows: In Section 2, the parameters' specification and their random combination techniques are described. Based on the random combination, the finite element modeling of the composite sheet is explained. In Section 3, the surrogate model construction and validation methods are given. Results and discussions are given in Section 4, and conclusions are given in Section 5.

Numerical Simulation
A carbon fiber reinforced polymer (CFRP) composite sheet with three plies, shown in Figure 1, was used to perform computer simulations. The dimensions and mechanical and thermal properties of this sheet are given in Table 1. The sheet was fixed at its left edge and free in the other edges. In the numerical simulation, the fact that during press forming, adjacent plies slip over each other was modeled. A uniform pressure load was applied on the surface of the sheet, and the middle ply was pulled out with different velocities.  To conduct the computer simulations, ABAQUS commercial software was used. Finite element modeling of the composite sheet was based on coupled temperature-displacement (C3D8T: an eight node thermally coupled brick, tri-linear displacement, and temperature) element type. The initial temperature (as taken from the autoclave) of the composite sheet was set to be 393.15 • K (120 • C). The fiber orientation angles of the two outer plies of the composite sheet were set to 0 • , relative to the ply pullout direction, throughout all the simulations, whereas that of the middle ply varied between 0 • and 90 • .
In the ABAQUS modeling of the inter-ply contact, the penalty scheme was used to impose frictional constraints on the tangential interaction of adjacent plies. This scheme allows the relative motion of plies while remaining stuck together. Moreover, the ply pullout speed is deliberately limited to be no more than 2 mm/s, as shown in Table 2. Under this very slow forming speed, the shear stress is not affected by inter-ply friction, and hence, the coefficient of friction was set to a constant value (0.1) in the forming simulation. The initial temperature for each computer simulation was set to 393.15 • K. During the forming process, for each sample data point, the temperature of the sheet (blank) drops down to its final value, which lies in the range given for the blank temperature parameter.

Parameters
The parameters were selected based on the reports of various researchers investigating the causes of composite forming defects. As mentioned in the Introduction, these parameters are: the relative fiber orientation angle of adjacent plies, forming temperature, pressure load, and ply pullout speed. The maximum and minimum limits of these parameters are given in Table 2. The selection of the minimum and maximum limits of the parameters was based on the thermo-mechanical properties of the composite material under consideration.

Design of Experiments
Each parameter was normalized between zero and one, then partitioned into a number of design points. The design points of the parameters were randomly combined using the Latin hypercube method. For surrogate model construction, forty-five sample data points were generated. Likewise, twenty-five validation data points were generated to validate the constructed model. Computer simulations were conducted for both sample and validation data points, and their corresponding maximum shear stresses were recorded.
The Latin hypercube method is one of the design of experiments methodologies that generates sample data points of the parameters. This Latin hypercube method partitions each of the design parameters and randomly combines them to generate sample data points. It is commonly preferred for computer based experiments [30][31][32]. Based on this methodology, the range of each parameter, given in Table 2, was divided into (45 for surrogate model training and 25 for model validation) non-overlapping intervals on the basis of equal probability.
Based on the probability density in an interval, a value from each interval was selected randomly. For instance, the 45 values of the fiber orientation angle obtained were paired randomly with the 45 values of the pressure load. These 45 pairs of the two parameters were then combined randomly with the 45 values of blank temperature to form 45 triplets, and then, with the values of ply pull-out speed, forty-five four-tuplets were generated. These 45 quadruplets were the same as the 45 four-dimensional input vectors whose component values were bounded in the ranges of the parameters shown in Table 2.

Surrogate Modeling
In many research works related to engineering disciplines, due to system complexity, the driving theoretical formulation for the system becomes difficult, or if derived, its numerical simulation could be computationally costly. Statistical approaches that take experimental data into account and model the theoretical formulation for the system are often applied and named "surrogate models". However, these models do not capture the exact functional form of the theoretical formulation, and approximations made by these models are viable, referring to the difficulty of the theoretical formulation or computational costs.
There are various types of surrogate models where polynomial regression, radial basis function, kriging, and neural networks are some of the commonly known ones. Yondo et al. [33] reported a detailed study on these surrogate models. There is no theoretical formulation for inter-ply shear stress as a function of the mentioned parameters. Therefore, in the work presented here, a function to be approximated is the inter-ply shear stress of the composite sheet, and the independent variables of the function are the mentioned forming parameters.

Surrogate Model
A surrogate model of inter-ply shear stress was developed. The developed model has a functional dependence on the aforementioned forming parameters. To approximate multidimensional functions, surrogate models such as the radial basis function are often used. The commonly known radial basis functions are given in Table 3.
In this study, multi-quadratic type radial basis functions (MQ-RBFs) are found to well approximate the functional form of inter-ply shear stress. Therefore, the MQ-RBFs are for the construction of the surrogate model. Table 3. Commonly known radial basis functions.

Radial Basis Function Mathematical Expression
Multi-quadratic

Radial Basis Function
The exponential convergence, dimension insensitivity, mesh, and geometry-free features of RBFs draw much attention in various research disciplines. These bases functions can be categorized into two: finitely and infinitely differentiable functions. Finitely differentiable RBFs depend only on Euclidean distance between two grid points, whereas infinitely differentiable RBFs depend on both the shape parameter and Euclidean distance between two grid points.
For the infinitely differentiable type of functions, shape parameter selection remains unsolved. An unrealistic value of this shape parameter affects the accuracy and stability of its RBF by either flattening (invoking the ill-possedness of the matrix) or sharpening (underestimating the effective arrangement of the control point). Radial basis functions have the following form, where φ(| x − x i | c) are user defined basis functions. The variables x, x i , and c imply interpolation points, control (center) points, and the shape parameter, respectively. The form of this basis function can be of any type as long as its positive definiteness is confirmed. Based on the type of RBF used, the shape parameter may or may not be required. λ i are weight parameters that can be obtained from the linearly independent system of equations.
The implemented MQ-RBFs have the form: where k stands for the number of independent design parameters and i represents the known data points, respectively. The approximating function and the inter-ply shear stress should have the same numerical values at least at the known data points. Based on this restriction, there are 45 independent equations where each equation is a linear combination of the non-linear MQ-RBFs. These equations can be written in matrix form as follows: The compact form of the above matrix is given as: where [Φ] −1 is the inverse of the square matrix. Based on the evaluated weight parameters, the 25 generated sample points, and the selected shape parameter, the functional form of inter-ply shear stress f (x) can be expressed as: wheref (x) is the surrogate model of f (x). This surrogate model can be used to predict the influence of a change in any of the independent variables (forming parameters). Therefore, this model is also known as the predictor. This surrogate model is developed based on a continuous and unidirectional carbon fiber reinforced thermoplastic composite sheet. Nevertheless, the model is applicable to all types of continuous and unidirectional fiber reinforced thermoplastic composite sheets.

Model Validation
Before the predictor is used for inter-ply shear stress analysis, its fidelity should be verified. Model accuracy is often validated by comparing output data obtained through computer simulation and the data obtained from the developed model. To validate the developed model, computer-based simulations were conducted for the 25 sample points, and their corresponding maximum shear stress (output) was obtained. The simulation results and the model outputs, using Equation (6), were compared to verify the fidelity of the model. The following commonly known statistical methods were used for the validation of the model.

•
Residual (RSD): the measure of the difference between the computer based simulation results and the model outputs. The mathematical form of RSD is, If the plot of the residual versus surrogate model shows a random distribution of points about the axis represented by the surrogate model, then this model is statistically considered as an appropriate model to be used for unknown function (inter-ply shear stress) approximation. The random distribution of points in Figure 2b shows that the model is valid for implementation. Furthermore, the linearity between the simulation results and model outputs, as shown in Figure 2a, tells us that the developed model is adequate to govern the evolution of the function to be approximated. • R-squared: the measure of how close the model versus simulation data are to the linear fit. The closeness of the data to the linear fit is shown in Figure 2a. The more the data are closer to the linear fit, the higher the value of R-squared is. The mathematical expression of R-squared is given as, If the data lie on the linear fit, then the developed model is said to be the exact functional form of the unknown function. The R-squared value of the developed model, as shown in Table 4, indicates that the data are close enough to the linear fit. However, the value of R-squared depends on the number of independent variables. R-squared increases with the increase in the number of independent variables. Therefore, the R-squared value is not always reliable. To resolve this reliability issue, often, the number of independent variables is incorporated into the R-squared, and the modified expression is termed the "adjusted R-squared". This is mathematically given as: where K represents the number of independent variables. • Root mean squared error (RMSE): measures the average model performance. The mathematical form of this measure is given as, The upper and lower limit of RMSE is zero and one, respectively. The value zero indicates the absence of residual variance. This means that the developed model is exactly equal to the unknown function. On the other hand, the value of one indicates that the developed model does not approximate the unknown function at all. RMSE is affected by the number of sample data points. To eliminate this effect, RMSE is normalized and hence named normalized root mean squared error (NRMSE), which is expressed as RMSE where f max and f min are the maximum and minimum values, respectively, of inter-ply shear stress (simulation results). Table 4. This error of 0.6 validates the adequacy of the model.

•
Standard deviation: measures the average amount of deviation of the predicted values from their mean. Its mathematical form is: The standard deviation is shown in Table 4.

Results and Discussions
This study aims to develop a surrogate model or predictor that approximates the functional form of inter-ply shear stress on four parameters: relative fiber orientation angle of adjacent plies, forming temperature, pressure load, and ply pullout speed. The fidelity of the developed model was verified through different statistical techniques.
The linear fit of predictor versus simulation, as well as residual versus predictor are shown in Figure 2a and Figure 2b, respectively. The closeness of the data points to the linear fit, in Figure 2a, shows that the model well predicts the functional dependency of inter-ply shear stress on the forming parameters. The horizontal line, in Figure 2b, shows the constancy of the average residual throughout the ranges of forming parameters. This in turn indicates that the accuracy of the developed model does not deteriorate over the ranges of forming parameters. Moreover, the random distribution (absence of pattern) of residual about the horizontal line shows that the surrogate model linearly fits the inter-ply shear stress. These statistical inferences indicate that the developed surrogate model is acceptable.
Other statistical methods shown in Table 4 also re-affirm the developed model's fidelity. The adjusted R-squared percentage value assures that the inter-ply shear stress data obtained through computer simulation is strongly related to that of the data obtained using the constructed model. The NRMSE of 0.6% indicates that there is moderate residual variance. The 1.71% average deviation of the values obtained via the developed model from the values obtained through computer simulations also confirms that the developed model well captures the functional form of inter-ply shear stress. After the developed model was confirmed to be valid to predict the functional dependency of the inter-ply shear stress, the model was used to analyze individual and paired effects of the forming parameters. The effects of individual parameter on inter-ply shear stress are shown in Figure 3. In these plots, all but one parameter are kept constant, as shown in Table 5, and predictor response is plotted against the free variable. Referring to Figure 3a,b, the overall reduction in inter-ply shear stress is observed as the relative fiber orientation angle of adjacent plies or forming temperature increases. The increase in pressure load hinders inter-ply slippage and hence increases inter-ply shear stress ( Figure 3c). As shown in Figure 3d, inter-ply shear stress increases with the increase in ply pullout speed up to, roughly, 0.8 mm/s and then reduces thereafter. The cause for the reduction of inter-ply shear stress could be an onset of ply slippage.  In Figure 4, surface plots of paired forming parameters versus inter-ply shear stress are shown. For each plot, two of the parameters were set to constant values, while the rest varied between their respective lower and upper bounds.
In Figure 4a, the values of the relative fiber orientation angle of adjacent plies and the pressure load were kept constant (refer to Table 5), while the values of blank temperature and ply pullout speed were set to vary. The plot shows that inter-ply shear stress decreases due to an increase in temperature and moderately increases as ply pullout speed increases. This supports the fact that inter-ply slippage is low at low blank temperature and hence increases inter-ply shear stress. In Figure 4b, the values of pressure load and blank temperature are set to their constant values, while the relative fiber orientation and ply pullout speed vary. The figure shows that inter-ply shear stress reduces as the relative fiber orientation angle of adjacent plies increases. Since the fiber orientation angles of the two outer plies are 0 • each, the maximum relative angle of the middle ply is 90 • . This affirms that the middle ply is easily deformable as ply pullout velocity is perpendicular to the fiber direction, and hence, the stress is low.
In Figure 4c, the pressure load and the ply pullout speed were set to their constant values, and the other parameters were set to vary. As relative fiber orientation angle and temperature increase, inter-ply shear stress decreases. The figure also depicts that the shear stress is more sensitive to changes in temperature than the relative fiber orientation of adjacent plies. The increase in pressure load and ply pullout speed, as shown in Figure 4d, increases inter-ply shear stress. High pressure load hinders inter-ply slippage and therefore increases inter-ply shear stress as shown in this figure.

Conclusions
To the best of authors' knowledge, there is no theoretically formulated mathematical model that governs the evolution of inter-ply shear stress during fiber reinforced composite sheet forming. Studies show that various parameters such as the blank temperature, pressure load, rate of forming, and relative fiber orientation angle of adjacent plies affect the formability of composite structures. In order to predict the influences of the mentioned forming parameters on inter-ply shear stress, a mathematical model is essential.
To develop the mathematical model of a function with multiple variables, the surrogate modeling approach is commonly preferred. Therefore, the surrogate model of inter-ply shear stress is developed based on numerical data obtained from computer simulation for various combinations of four forming parameters. Through various statistical methods, the developed surrogate model is then proven to well approximate the functional dependency of inter-ply shear stress on these parameters.