Approximation of Non-Linear Stress–Strain Curve for GFRP Tensile Specimens by Inverse Method

Studying the characteristics of materials through a finite element analysis (FEA) has various benefits; hence, many studies have been conducted to improve the reliability of the analysis results. In general, the mechanical properties used in FEA for metals and metal composites are stress–strain data obtained through tensile tests, which are used for modeling from a macroscopic perspective. While many studies have been conducted on metal materials, there are limited studies on the analysis of polymer composite materials produced through injection and special processing. In this study, existing inverse methods were applied, and an FEA was conducted to reproduce the axial displacement of the tensile specimens comprising glass fiber-reinforced polymer (GFRP); further, errors were examined by comparing the test and analysis results. To reduce such errors, the experiment and the FEA results were analyzed through parameter optimization based on various empirical formulas. The accuracy of various inverse methods were examined and an inverse method suitable for GFRP was proposed.


Introduction
Many studies have been conducted to reduce the weights of automotive [1][2][3][4] and aircraft parts [5,6]. A few of these parts are made of metals; replacing them with polymer plastics satisfies the demand for weight reduction [7,8]. However, polymer plastics exhibit lower strength compared with metal materials [9,10]. To this end, many studies have been conducted to improve the strength of plastics by adding glass or carbon fibers [11,12]. The plastics with such fibers vary the pattern of the stress-strain (S-S) curve, a macroscopic property, and also make it difficult to predict microscopic mechanical properties [13,14]. Unal [15] proposed a method for estimating the S-S curve of a specific or composite material with a large fluctuation in material properties during processing, through fuzzy theory. As there is a considerable difference between the S-S curve derived through fuzzy theory and the behavior of an actual specimen, it is difficult to use the derived S-S curve as the material property for a finite element analysis (FEA).
Recently, studies on models for estimating the S-S curve using the response surface method of design of experiment (DOE), which are referred to as inverse methods, have been conducted [16][17][18]. For inverse methods, basic theories are provided by the American Society for Testing and Materials (ASTM) standards; further, studies to improve the accuracy of property data through mathematical modeling that uses power and exponential functions have been proposed [12][13][14][15][16][17][18][19][20][21][22].
However, as mathematical models may vary depending on the materials, inverse methods have been mostly applied to metals and alloy models. Therefore, it is necessary to select and optimize inverse methods to secure the macroscopic properties of composite materials.
In this study, a method for estimating the macroscopic properties of glass fiber-reinforced polymer (GFRP) tensile specimens, which are commonly used as the lightweight materials of automotive parts, using an inverse method through the finite element method (FEM) was proposed. GFRP specimens were fabricated, and their mechanical properties for tensile deformation were analyzed through tests. Based on the results, various inverse methods were applied for modeling, and the analysis and experiment results were compared. Moreover, parameter optimization was conducted using DOE to select an inverse method most suitable for GFRP among various models. In this process, the maximum number of iterations for optimization was limited to the same value. Moreover, whether the actual properties of GFRP could be estimated was verified through advanced inverse methods-by combining the function with the lowest error.

Experimental Setup (ASTM D638-14)
In this experiment, SUPRAN ® PP 1320 (wt 20% GF), 1330 (wt 30% GF), and 1340 (wt 40% GF) of SAMBARK LFT were used. Each material is based on polypropylene, and 6~25 mm-long glass fiber accounts for 20%, 30%, and 40% of the total weight of the material. Figure 1a shows the conceptual diagram of the injection process, and Figure 1b shows the test specimen produced by the injection process. The process conditions were set to the injection temperature of 220 • C, pressure of 50 MPa, mold temperature of 4 • C, and cooling time of 15 s, respectively. ASTM D638-14 [23] was applied for the tensile test method. The specimens' geometry was the dog-bone type, as shown in Figure 1c, and they were fabricated by injection molding, as shown in Figure 1b. To observe the various tendencies of the inverse analysis for GFRP, specimens were fabricated using the properties of GFRP wt 20%, 30%, and 40% GF. Three tensile specimens were fabricated for each GFRP property.
Tensile tests on the specimens were conducted using the MTS 810 material test system. Figure 2 shows the schematic of the tensile test setup. MTS 634 12F was used as the strain gauge equipment for measuring gauge length (G), as shown in Figure 1c. Tensile tests were conducted considering a cross-head speed of 5 mm/min, in accordance with ASTM D638-14. The variations of the load and gauge length were obtained as data and reflected to the engineering stress and strain. Figure 3 shows the properties obtained through tensile tests. For each property, three tests were conducted, and the data corresponding to the median values were selected as representative values. Force (F a ) derived from the tensile tester and displacement (∆G) derived from the strain gauge were used to express the engineering strain-stress data as σ e and ε e as shown below, along with the true stress presented by ASTM D638-14: Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 18 applied to metals and alloy models. Therefore, it is necessary to select and optimize inverse methods to secure the macroscopic properties of composite materials. In this study, a method for estimating the macroscopic properties of glass fiber-reinforced polymer (GFRP) tensile specimens, which are commonly used as the lightweight materials of automotive parts, using an inverse method through the finite element method (FEM) was proposed. GFRP specimens were fabricated, and their mechanical properties for tensile deformation were analyzed through tests. Based on the results, various inverse methods were applied for modeling, and the analysis and experiment results were compared. Moreover, parameter optimization was conducted using DOE to select an inverse method most suitable for GFRP among various models. In this process, the maximum number of iterations for optimization was limited to the same value. Moreover, whether the actual properties of GFRP could be estimated was verified through advanced inverse methods-by combining the function with the lowest error.

Experimental Setup (ASTM D638-14)
In this experiment, SUPRAN ® PP 1320 (wt 20% GF), 1330 (wt 30% GF), and 1340 (wt 40% GF) of SAMBARK LFT were used. Each material is based on polypropylene, and 6~25 mm-long glass fiber accounts for 20%, 30%, and 40% of the total weight of the material. Figure 1a shows the conceptual diagram of the injection process, and Figure 1b shows the test specimen produced by the injection process. The process conditions were set to the injection temperature of 220 °C, pressure of 50 MPa, mold temperature of 4 °C, and cooling time of 15 s, respectively. ASTM D638-14 [23] was applied for the tensile test method. The specimens' geometry was the dog-bone type, as shown in Figure 1c, and they were fabricated by injection molding, as shown in Figure 1b. To observe the various tendencies of the inverse analysis for GFRP, specimens were fabricated using the properties of GFRP wt 20%, 30%, and 40% GF. Three tensile specimens were fabricated for each GFRP property.
Tensile tests on the specimens were conducted using the MTS 810 material test system. Figure 2 shows the schematic of the tensile test setup. MTS 634 12F was used as the strain gauge equipment for measuring gauge length (G), as shown in Figure 1c. Tensile tests were conducted considering a cross-head speed of 5 mm/min, in accordance with ASTM D638-14. The variations of the load and gauge length were obtained as data and reflected to the engineering stress and strain. Figure 3 shows the properties obtained through tensile tests. For each property, three tests were conducted, and the data corresponding to the median values were selected as representative values. Force (F a ) derived from the tensile tester and displacement (∆G) derived from the strain gauge were used to express the engineering strain-stress data as σ e and ε e as shown below, along with the true stress presented by ASTM D638-14: (2) (3) Table 1 lists the Young's modulus (E), Poisson's ratio (ν), yield stress (σ y ), and tensile strength (σ ) derived through tensile tests.    Table 1 lists the Young's modulus (E), Poisson's ratio (ν), yield stress (σ y ), and tensile strength (σ UTS ) derived through tensile tests.

Finite Element Analysis Setup
Figure 4a-e shows the finite element models of the tensile specimen. In the models, the width (W, W 0 ) of the tensile specimen was divided into two, four, six, eight, and ten elements, while the aspect ratio of the elements ranged from 0.5 to 1.5. As the geometry of the specimen was not complex, two-dimensional (2D) triangle elements and three-dimensional (3D) tetra elements were not considered.
As for the material properties used for the conformity evaluation of the elements, the experimental parameter (σ e ) of GFRP wt 20% is commonly applied. The finite element models were analyzed under the conditions shown in Figure 4f based on ASTM D638-14. The left section of the specimen model was set as a fixed section, while the right section was moved at a speed of 5 mm/min. When the FEA was conducted, the reactive force was obtained from the left node set of the specimen model, and the displacement of the gauge length was obtained using the specimen information.

Finite Element Analysis Setup
Figure 4a-e shows the finite element models of the tensile specimen. In the models, the width (W, W 0 ) of the tensile specimen was divided into two, four, six, eight, and ten elements, while the aspect ratio of the elements ranged from 0.5 to 1.5. As the geometry of the specimen was not complex, two-dimensional (2D) triangle elements and three-dimensional (3D) tetra elements were not considered.
As for the material properties used for the conformity evaluation of the elements, the experimental parameter (σ e ) of GFRP wt 20% is commonly applied. The finite element models were analyzed under the conditions shown in Figure 4f based on ASTM D638-14. The left section of the specimen model was set as a fixed section, while the right section was moved at a speed of 5 mm/min. When the FEA was conducted, the reactive force was obtained from the left node set of the specimen model, and the displacement of the gauge length was obtained using the specimen information. The implicit solver and explicit solver of LS-Dyna R910 were used as FEA tools. The analysis results are shown in Figure 5 according to the division of elements, solver type, and dimensions. Figure 5a shows that the division of elements and solver type did not significantly affect the accuracy of the analysis. However, they significantly affected the analysis time, as shown in Figure 5b. The implicit solver and explicit solver of LS-Dyna R910 were used as FEA tools. The analysis results are shown in Figure 5 according to the division of elements, solver type, and dimensions. Figure 5a shows that the division of elements and solver type did not significantly affect the accuracy of the analysis. However, they significantly affected the analysis time, as shown in Figure 5b In general, models with implicit solvers, 3D elements, and a large number of divisions are selected to improve the accuracy of the analysis results. However, when estimating the mechanical properties, it was confirmed that selecting a model with explicit solvers, 2D elements, and a small number of divisions did not result in large differences. Therefore, in this study, models with explicit solvers, 2D elements, and a small number of divisions were used to reduce the analysis time.

Inverse Methods
Inverse methods are useful when the change in the thickness of the necking zone cannot be accurately measured after the yield stress is reached. They can be used as correction measures for the stress value when the engineering S-S curve obtained experimentally is different from the engineering S-S curve obtained using the FEA. In general, models with implicit solvers, 3D elements, and a large number of divisions are selected to improve the accuracy of the analysis results. However, when estimating the mechanical properties, it was confirmed that selecting a model with explicit solvers, 2D elements, and a small number of divisions did not result in large differences. Therefore, in this study, models with explicit solvers, 2D elements, and a small number of divisions were used to reduce the analysis time.

Inverse Methods
Inverse methods are useful when the change in the thickness of the necking zone cannot be accurately measured after the yield stress is reached. They can be used as correction measures for the stress value when the engineering S-S curve obtained experimentally is different from the engineering S-S curve obtained using the FEA.
Inverse methods to obtain the true stress are largely divided into power law models and exponential law models. The models of Hollomon [24], Ludwik [25], Swift [26], and Gosh [27] shown in Equations (4)- (7) can be considered as representative power law models. The models of Voce [28] and Hockett-Sherby [29] shown in Equations (8) and (9) can be considered as representative exponential law models. In each model, c is a material parameter, and its initial value is determined using the logarithmic curve method given in ASTM E646-16 [30].
In this theoretical model, parameters of the test environments such as the ambient temperature and test speed can be included. However, in this study, the parameters of the test environments were not considered, as the aim of the study was to propose an inverse method that could approximate the plastic deformation behavior of GFRP. Figure 6 shows the true S-S curves derived using Equations (4)- (9). In the figure, test data for the curves were included for a comparison between the true stress and the engineering stress.
It was found that the true S-S data obtained using the ASTM D638-14, Hollomon, Ludwik, and Swift equations were not significantly different from the experimental data. In contrast, the Gosh model and the exponential law model exhibited significant differences from the experimental data. The physically reliable and accurate data of the engineering S-S curve and all true stress curves were estimated values.
The true S-S curve data were used as the mechanical properties of the FEA software, and were entered as the point data of 100 divisions based on the total strain. The engineering S-S curve data shown in Figure 7 were derived using the simulation shown in Figure 4f.
To quantitatively demonstrate the expression accuracy of the analysis results for an actual physical phenomenon, the engineering S-S curve obtained experimentally was compared with the engineering S-S curve obtained from the analysis, as shown in Equation (10). Through this process, an inverse analysis model that could accurately express actual physical phenomena was found.
The difference in the engineering S-S curves obtained using the experimental and analysis data was expressed as e q , as shown in Equation (10), and the average of the absolute values of such data was expressed as e q , as shown in Equation (11). The average of the errors generated from all of the properties for each inverse method was expressed as e r , as shown in Equation (12).    Figure 7 shows the results of the FEA that was conducted using the estimation equations. For the engineering S-S curves, the experimental data were significantly different from the analysis data, and the c value of ASTM E646-16 resulted in inaccurate results. This can also be observed in Figure  8, which shows the graphs obtained using Equation (10) for a comparison. The arithmetic average values for e q are listed in Table 2. Among the power law models, the Gosh model exhibited the most accurate approximations of the deformation when compared to the actual values, and among the exponential law models, the Voce model exhibited the best approximation.  Figure 7 shows the results of the FEA that was conducted using the estimation equations. For the engineering S-S curves, the experimental data were significantly different from the analysis data, and the c value of ASTM E646-16 resulted in inaccurate results. This can also be observed in Figure 8, which shows the graphs obtained using Equation (10) for a comparison. The arithmetic average values for e q are listed in Table 2. Among the power law models, the Gosh model exhibited the most accurate approximations of the deformation when compared to the actual values, and among the exponential law models, the Voce model exhibited the best approximation.

Optimization of Inverse Methods
To approximate the experiment and analysis results for the engineering S-S curve, the value of the material parameter (c) of the inverse method must be adjusted. To this end, optimization was performed using the minimum error between the experimental and the analysis results as the objective function.
The optimization formulation process is shown in Equations (13)- (18). To improve the accuracy of the engineering S-S curve obtained using the FEA, an objective function that minimizes e q was constructed, and the material parameter was adjusted in the 30% range from the initial value used given in ASTM E646-16. Find: x ∈ R 3 (13) Minimize: Subject to: x k_min ≤ x k ≤ x k_max ; k = 1, 2, 3 . . . K Parameter Optimization D-optimal DOE was applied to minimize the number of iterations in the optimization process [31]. Here, the process of deriving the error range between the experimental and FEA results and adjusting the parameters was repeated. The same optimization method and number of iterations were applied to all inverse methods. Figure 9 shows the analysis results. As the engineering stress and true stress models based on ASTM D638-14 were excluded from the optimization targets, the errors indicated by Figure 7 could be observed. When the optimized inverse methods were used in the FEA, the results were very similar to the engineering S-S curve obtained from the experiment. The error of the optimized inverse method (e q ) can be seen in Figure 10. The error patterns of the power law and that of the exponential law models were reversed, and it can be seen that the Hockett-Sherby formula was very accurately optimized.
The average of cumulative errors (e q ) is listed in Table 3. Considerable errors were observed in the three data points excluded from the optimization process. Among the power law models, the Swift model obtained the most accurate results, followed by the Ludwik, Gosh, and Hollomon models. Among the exponential law models, the Voce model obtained the most accurate results, followed by the Hockett-Sherby model.    Figure 10. Error values of the optimized inverse method models: (a) optimized inverse models of GFRP wt 20% GF; (b) optimized inverse models of GFRP wt 30% GF; (c) optimized inverse models of GFRP wt 40% GF.

Advanced Inverse Method
It is known that power law models are suitable for metals with the BCC structure, such as steel, and that exponential law models are suitable for metals with the FCC structure, such as aluminum and copper. In addition, Kim-Tuan [32] combined the power law (Swift) model and the exponential law (Voce) model for the efficient management of various metal properties, as shown in the following equations. They proposed advanced inverse methods to adjust the proportion of each model using material parameters.
These advanced inverse methods are compatible with various properties, and can correct the errors of the power and the exponential laws. Figure 10 shows that the error (e q ) values of the power laws and exponential laws optimized in this study had opposite patterns. Therefore, if the accurately derived power and exponential laws are selected and combined into one equation, better accuracy can be obtained. Through the optimization results, the two models with the lowest average property error (e r ) values were selected from the power and exponential laws.
The Swift and Ludwik models from the power laws and the Voce and Hockett-Sherby models from the exponential laws were combined to use the advanced inverse methods of the following equations. In each advanced inverse method, "d" is the material parameter that adjusts the proportions of the power and exponential laws. Figure 11 shows the e q values, including the optimization results of the advanced inverse methods. Items 8-12 in Table 3 are the optimized results of the advanced inverse methods. It can be seen that the advanced inverse methods can estimate the actual behavior of the specimen more accurately. By comparing the optimized e r , it is seen that the Swift-Voce model exhibited the most approximate estimate of 0.18 MPa. When the levels of factors increased in the DOE, the order of the data analysis also increased; hence, the non-linear response must be considered. Therefore, the number of factors must be reduced to increase the precision and decrease the cost [33]. In many studies, the number of factors was reduced, and the DOE as well as the optimization was performed for the simplified models. In addition, their performances were verified [34]. As the advanced inverse methods have more factors than the normal inverse methods, it is difficult to obtain high accuracy under the same DOE and with the same number of optimization iterations. The Hockett-Sherby model of the exponential law exhibited the most accurate value.

Simulation and Optimization
However, the advanced inverse models exhibited higher accuracy compared with the inverse models, even though they had two to three times more factors (Figure 11). This indicates that plastic and composite materials are significantly different from metals, which can be expressed by one power or exponential function.
under the same DOE and with the same number of optimization iterations. The Hockett-Sherby model of the exponential law exhibited the most accurate value.
However, the advanced inverse models exhibited higher accuracy compared with the inverse models, even though they had two to three times more factors (Figure 11). This indicates that plastic and composite materials are significantly different from metals, which can be expressed by one power or exponential function. Figure 11. Error values of the optimized inverse method. Figure 11. Error values of the optimized inverse method.

Conclusions
In this study, the experiment results were compared with the results of the analysis, which was conducted by applying the proposed empirical equations to examine the S-S curve of GFRP tensile specimens through a FEA. In the study, inverse models were defined using the coefficients based on ASTM E646-16, and the results of the FEA were found to be significantly different from the experiment values. Among the models by the logarithmic curve, the Voce model exhibited the smallest difference of less than 1.5 ± 0.5 [MPa]. To improve the accuracy of the analysis results, the factors of the inverse methods were optimized, and curve fitting was performed. As a result, a relatively low error rate was observed, compared with the existing models. After optimization, the Hockett-Sherby model exhibited the lowest error rate. These results indicate that the application of the Hockett-Sherby model can improve the accuracy of the analysis results when a FEM analysis is conducted using the macroscopic mechanical properties of GFRP.   Difference in the experiment and analysis for the engineering stress e q Average of the differences between the experiment and the analysis e r Average of the differences between the experiment and the analysis for the entire material properties c Material parameter of inverse model d Material parameter of advanced model K Maximum number of material parameters by inverse method k Value that expresses the order of a material parameter