A Phenomenological Primary–Secondary–Tertiary Creep Model for Polymer-Bonded Composite Materials

This study develops a unified phenomenological creep model for polymer-bonded composite materials, allowing for predicting the creep behavior in the three creep stages, namely the primary, the secondary, and the tertiary stages under sustained compressive stresses. Creep testing is performed using material specimens under several conditions with a temperature range of 20 °C–50 °C and a compressive stress range of 15 MPa–25 MPa. The testing data reveal that the strain rate–time response exhibits the transient, steady, and unstable stages under each of the testing conditions. A rational function-based creep rate equation is proposed to describe the full creep behavior under each of the testing conditions. By further correlating the resulting model parameters with temperature and stress and developing a Larson–Miller parameter-based rupture time prediction model, a unified phenomenological model is established. An independent validation dataset and third-party testing data are used to verify the effectiveness and accuracy of the proposed model. The performance of the proposed model is compared with that of an existing reference model. The verification and comparison results show that the model can describe all the three stages of the creep process, and the proposed model outperforms the reference model by yielding 28.5% smaller root mean squared errors on average.


Introduction
Creep is a time-dependent progressive inelastic deformation behavior. It can cause the relaxation of stress and irreversible deformation, thus leading to functional failures when the part is intended to maintain the required stress and shape [1]. Polymer-bonded composites materials (PBMs) are increasingly used in engineering components due to its high strength and lightweight [2][3][4]. For PBMs the creep can take place at a relatively low temperature [5]; therefore, a reliable prediction of the creep behavior is critical to ensure the safety and durability of PBMs under sustained loads.
Several studies have shown that the creep behavior of PBMs depends on many factors such as matrix content [6][7][8][9], particle size [10][11][12], stress [13][14][15], temperature [16,17], and humidity [18,19]. The importance of long-term properties for polymer composites are highlighted in recent studies [20,21]. Experimental results reveal that the entire process of creep deformation can be divided into three stages [22,23], namely the primary (transient) creep, the secondary (stationary) creep, and the tertiary (unstable) creep. The three stages of the creep process are illustrated in Figure 1a with the corresponding creep strain rate shown in Figure 1b. The secondary (stationary) creep is considered to be the dominant creep for many applications. In this stage, the equilibrium between the softening and hardening of the material is assumed, leading to a stable strain rate [24][25][26]. Prior to the stationary stage, a short transient period of primary creep is required to reach such an approximate equilibrium between the softening and hardening processes. The final part of the creep process is the tertiary stage where the strain rate increases rapidly until rupture [27,28]. The progressive damage such as the formation and growth of voids on grain boundaries are considered to be the main contributors of the rapid growth of the strain rate. Traditionally the three stages are modeled separately. For primary creep, the creep rate equations are in general given by [29] ε = dε dt = aσ n · t m Time hardening aσ n · ε m Strain hardening , where a, n, and m are fitting parameters. For the secondary creep, the equivalent creep rate equations are given by Equation (2) [30][31][32][33][34][35][36], where a, b, a 1 , a 2 , σ 0 , n, n 1 , and n 2 are fitting parameters, and σ is the effective stress, e.g., the equivalent stress. To accommodate the temperature effect, the Arrhenius law is usually incorporated into Equations (1) and (2). Existing constitutive modeling methods for the creep behavior under different testing conditions can be classified into three categories, the general stress-strain-time modeling, the rheological modeling, and the empirical modeling [37]. The generalized stress-straintime modeling is used to describe the viscous effect and rate-independent behavior under general loading conditions, including the creep model based on the elastic-viscoelastic correspondence principle and a material stiffness equation [38], the Wiechert model with damage evolution law and time-temperature shift factor for creep response prediction under different temperature [39], the micromechanics model employing the correspondence principle in viscoelasticity [40], and the viscoelasticity-viscoplasticity temperature-dependent model that includes anisotropic damage evolution [41]. The rheological modeling is centered on the idea of using basic components such as springs, dashpots, and sliders to describe the creep strain behavior with time. The creep model of PBM can be formulated by combining various basic components with the viscoelastic principle. For example, the Burgers model based on time-temperature equivalence principle [28,42], the fractional Poynting-Thomson model [43], and the viscoelastic response-based models [44,45]. The empirical modeling is to establish the correlations between the creep rate and other independent variables such as temperature and stress using testing data. Examples of this type of model include, but are not limited to, the stress-dependent phenomenological viscoelastic-plastic model [46], the modified power law model with temperature-and stress-dependent correction factors [47], and the improved Findley-Khosla model with the Schapery's integral [48]. For empirical modeling, the detailed mechanisms of the creep behavior are not fully explained. However, when the boundary conditions are consistent with the experiment, it can provide solutions for practical engineering problems [49].
The models for predicting the full-stage creep behavior of PBMs are limited. Krankel et al. developed a rheological model by introducing the time variable into the Burgers model for creep behavior prediction of bonded anchors [50]. Sudduth et al. developed a polynomial strain rate model which enables the capture of the behavior of the full creep curve [51]. However, the above two models cannot directly incorporate the dependence of model parameters on temperature and holding stress, making it difficult to predict the creep strain under more general conditions without testing data. The purpose of this study is to develop a unified phenomenological creep model for polymer-bonded composite materials, allowing for prediction of the creep behavior and life under more general conditions without testing data. To achieve that, the temperature-and stress-dependent effect is incorporated in the proposed rational function-based full creep strain rate equation, and the creep rupture time model based on the Larson-Miller parameter is developed. The former creep rate equation differs from the existing full-stage creep rate models in format and uses less (three) fitting parameters. The latter rupture time model under general conditions without testing data is rarely seen for PBMs.
The remainder of the study is organized as follows. First, the experimental work and creep data of a typical PBM are presented. Uni-axial compression creep testing in the temperature range of 20 • C-50 • C and the holding stress range of 15 MPa-25 MPa is performed to obtain the full-stage creep data. Next, a rational function-based phenomenological model is proposed to describe the creep behavior in the entire primary-secondary-tertiary process. By correlating the fitting parameters with temperature and stress, a general temperatureand stress-dependent full-stage creep model is formulated. The creep rupture time model is developed using the Larson-Miller parameter. Following that, the model is validated using independent testing data and third-party testing data of another type of PBM. The performance of the proposed model is further compared with an existing reference model. Finally, conclusions are drawn based on current results. Figure 2 presents the procedure of the overall methodology development of this study. In the experimental part, a total number of six specimens are prepared for uni-axial compression creep testing at temperatures ranging from 20 • C to 50 • C with a holding stress of 15 MPa-25 MPa. The total strain vs. time data are acquired for each of the specimens, and the strain rate data are extracted using a central difference scheme for strain rate model development.

Specimens and Experimental Setup
The PBM used in this study consists of 94 wt.% barium sulfate grains as filler material and 6 wt.% fluororubber as matrix material. The size of the filler is in the range of 0.5 mm and 3 mm. Cylinder specimens with dimensions shown in Figure 3 are prepared according to the standard [52].
Thirdparty data

Reference model
Temperature and stress dependent a(T, Strain rate data extraction   Prior to creep testing the specimens are examined using cone-beam computed tomography (CT) to ensure no initial damage exists in the materials. The inspection process is illustrated in in Figure 4, and a typical CT image is given where the light-colored filler particles are barium sulfate, and the rest dark area is the fluororubber binder material. A total number of six specimens are prepared, and testing conditions for the specimens are shown in Table 1. The uni-axial compression creep testing is performed using a universal testing machine with an environmental chamber. The environmental chamber allows for keeping the temperature at a prescribed value during the creep testing. The compressive stress is applied to the specimen at a strain rate of 0.5 mm/min until the prescribed stress is reached. After that the applied stress is sustained until rupture.

Creep Testing Results
The strain vs. time data of the six specimens from the initial state to rupture are acquired, and results are presented in Figure 5. It can be observed from Figure 5 that the creep behavior of the PBM exhibits three distinct stages. The initial rapidly increasing of the strain is mainly due to elastic deformation, plastic deformation, and work hardening.
After that, the strain curve remains a relatively constant slope for a significant amount of time. Following that is another rapidly increasing of the strain leading to the final rupture. In addition, a higher holding stress can greatly reduce the rupture life as shown in Figure

Strain Rate Extraction
A central difference scheme given as Equations (3) and (4) is used to extract the strain rate from the strain vs. time data.
where dε/dt is creep strain rate and the t is the time variable. The subscript i is the data point index and the total number of data points is n. The strain vs. time testing data are processed using Equations (3) and (4) to obtain the strain rate data. The extracted creep strain rate data of the five specimens used for model development are shown in Figure 6. The results presented in Figure 6 show that the primary and tertiary stages are much shorter than the secondary stages, and creep rate in the entire creep process varies by up to several order of magnitudes. Therefore, direct modeling of the creep strain rate data in linear scale may yield undesired rounding errors due to such large differences in data. On the other hand, the strain variation is monotonic which ensures that the strain rate is positive, allowing for logarithm transformation. In this study, the log-transformed strain rate data are used for model development.

Creep Model Development
The log-transformed creep strain rate data under each of the testing conditions resemble a bathtub shape; therefore, a rational function-based equation is proposed to model the log-transformed strain rate data. The resulting fitting parameters of the model are subsequently correlate with temperature and stress using a response surface model to incorporate the effects of temperature and stress.

Creep Strain Rate Model
The following rational function-based equation is proposed to describe the logtransformed strain rate data.
where dε/dt and t are defined as before, and a, b, and c are the fitting parameters. The parameter a loosely measures the magnitude of the strain rate in the secondary stage, and the terms b and c are related to the transitions from the primary to the secondary and the secondary to the tertiary, respectively. The term t ini is the initial time of the creep process, and t rup is the creep rupture time.
The fitting parameters using data associated with the five specimens are obtained using the regular nonlinear least square estimator, and the resulting model parameters are presented in Table 2. With the fitting parameters, the mean curves of the proposed equation are computed and shown in Figure 7. It can be observed that the proposed creep strain rate equation Equation (5) can reliably capture the three stages of the actual strain rate data. Table 2. Results of model fitting parameters (a, b, c) using Equation (5). The initial and rupture times are directly obtained from the raw testing data. The numbers in the first column corresponds to the test condition in Table 1.

No.
T Exponentiation of Equation (5) to recover the strain rate in linear scale aṡ The creep strain at a given time t can be obtained by time integration of Equation (6) from the initial time t ini to t as where ε 0 is the transient elastic-plastic strain which can be set as a prescribed value in testing or estimated using stress-strain constitutive models [53]. It is noted that Equation (7) can be resolved using numerical integrators such as RK45 and its variants with the initial value of ε 0 . Using Equation (7) and model parameters in Table 2, the strain vs. time results are obtained for the five specimens. Model results and the actual raw creep strain data are presented in Figure 8, where a close agreement between the two is observed.

Temperature and Stress Dependence Model
Based on the resulting model parameters under each of the conditions shown in Table 2, a first-order response surface model is employed to correlate the parameters with temperature and stress. The response surface model can be expressed as Equation (8).
where the α i , β i , and γ i (i = 1, 2, 3) are fitting coefficients, T is the temperature, and σ is the stress. Using the data in Table 2, the fitting coefficients α i , β i , and γ i , i = 1, 2, 3 are obtained using the regular least square estimator as The prediction results of the response surface model with the fitting coefficients given in Equation (9) are evaluated and presented in Figure 9.  Combining Equation (5) and Equation (8) with the parameters in Equation (9), the temperature-and stress-dependent creep strain rate model can be expressed aṡ The histogram of the model residuals is presented in Figure 10. The standard deviation of the residuals is estimated as 0.6382. --

Creep Rupture Time Model
It is noted that the rupture time variable t rup is required in the developed strain rate model Equation (10). The rupture time is defined as the time duration between the time when the part is loaded with the sustained stress and the time of the final fracture. In Equation (10), the rupture time can alter the tail region behavior of the strain curve. For conditions without tested specimens the corresponding rupture time is unknown. Consequently, the existing data on specimens tested under uni-axial compression are not sufficient for more general loading conditions. To predict the strain rate response under other conditions, it is necessary to establish rupture time prediction model to obtain t rup for a given stress and a temperature.
In this study, the Larson-Miller parameter (LMP) [54] is adopted to develop the rupture time prediction model of PBMs. LMP is an equation to calculate the creep rupture time at different temperatures under a given stress. The basic form of the LMP can be expressed as where T is the temperature in Celsius, t rup is the creep rupture time, and C and m are material constants, respectively [54]. To further introduce the stress variable into the rupture time prediction, a linear relationship between the stress and the LMP is proposed as where p 0 and p 1 are fitting parameters. Using the rupture time, temperature, and stress data in Table 2, the optimal parameters of (C, m) are identified as (50.998, 0.549) using the nonlinear least square estimator, and parameters p 0 and p 1 in Equation (12)  The rupture time under for a given combination of temperature and stress is obtained by substituting Equation (11) into Equation (12) as Incorporating Equation (13) into Equation (10) to obtain the final creep strain rate model as,ε By further substituting Equation (14) into Equation (7) to have the final creep strain model The model prediction results using Equation (15) and the actual creep strain testing data are compared in Figure 12. To quantify the performance of the model, the root mean squared error (RMSE) defined in the following equation is employed.
where y i is the actual value,ŷ i is the prediction value, and i = 1, . . . , N represents the index of a total number of N data points. RMSEs of the testing data on the five specimens used for model development are calculated using Equation (14) and presented in Table 3. The maximum RMSE is 0.0027 for the testing data obtained under the condition of T = 30 • C and σ = 20 MPa.

Model Validations and Comparisons
An independent dataset, Specimen 5 in Table 1, is used to validate the performance of the model. Moreover, the third-party testing data reported in Ref. [55] are used to validate the effectiveness of the model for other PBMs. In addition, the proposed model is compared with an existing reference model to demonstrate its performance, and the performances of the two models in terms of RMSE are compared and quantified.

Model Validation
Testing data of the validation specimen (No. 5 in Table 1) are used for validation. Using Equations (14) and (15) with the corresponding temperature and stress of the specimen, the creep strain results are obtained. Comparison between the model prediction results and the actual data are presented in Figure 13. It can be seen that the model can effectively capture the three stages of the creep process. To further investigate the generality of the proposed model, the third-party testing data on high-density polyethylene (HDPE) are employed [55]. The testing data consist of creep strain results of a total number of seven specimens. Model parameters of Equation (5) are obtained using data associated with the seven specimens. The response surface model coefficients are subsequently obtained using Equation (10). The resulting temperature-and stress-dependent model parameters for the third-party testing data are Using Equation (15) and the above model parameters, the creep strain prediction results are obtained and presented in Figure 14. A general close agreement between the model prediction results and the actual testing data can be observed, indicating the proposed model can be effectively used for creep strain prediction for other PBMs with similar strain behaviors. RMSEs of the model prediction results are evaluated and shown in Table 4 where the maximum value is 0.1196 under the condition of 53 • C and 8.8 MPa. Table 4. RMSEs of the model prediction for the creep strain data reported in Ref. [55].

Model Comparisons
To further demonstrate the performance of the proposed model, the model is compared with a reference model reported in Ref. [51]. The reference model in Equation (18) can also describe the primary-secondary-tertiary creep behavior. dε dt = εP 1 t · 1 + P 2 (P 3 ε) + P 4 (P 3 ε) 2 1 + P 1 + (2 + P 1 )P 2 P 3 ε + (3 + P 1 )P 4 (P 3 ε) 2 , where P 1 , P 2 , P 3 , P 4 are fitting parameters and other variables are defined as before. The same data in Table 1 for modeling are used to obtain the required model parameters of Equation (18). The prediction results of creep strain are evaluated using the proposed model and the reference model, and are compared with the actual testing data in Figure 15. In general, the two models both yield satisfactory fitting results.
RMSEs of the two models under each of the conditions are evaluated and compared in Figure 16, where the proposed model yields smaller RMSEs in all data sets except for the case of (30 • C, 20 MPa). The sum of the RMSEs associated with the proposed model and the reference model are 0.73 × 10 −3 and 1.02 × 10 −3 , respectively. The proposed model reduces the overall RMSE by about 28.5%.
In addition, the proposed model predicts the tertiary creep stage more reliably than the reference model. Fluctuations of the resulting creep strains produced by the reference model can be observed when the creep strain approaches the rupture life, as shown in Figure 16. For the strain rate equation of Equation (18), there exists a critical strain larger than which a negative strain rate can be produced by the equation. The negative strain rate reduces the creep strain to a value lower than the critical strain, leading to a positive strain rate again. This alternating nature of the polynomial function causes the strain rate varies between the negative and positive values. Consequently, the resulting strain fluctuates around the critical strain.

Conclusions
A unified phenomenological creep model was developed for polymer-bonded composite materials, allowing for predicting the creep behavior in the entire primary, secondary, and tertiary stages. A total number of six specimens made of a typical polymer-bonded composite material were prepared. The uni-axial compression creep testing with the holding stress in the range of 15 MPa-25 MPa at a temperature ranging from 20 • C to 50 • C was performed to acquired creep strain data. A rational function-based equation was proposed to describe creep strain rate in the entire creep process. The model parameters were identified using testing data. The temperature-and stress-dependent effect was incorporated into the strain rate model using a first-order response surface model. The creep rupture time model based on the Larson-Miller parameter is established, allowing for predicting the creep strain under more general conditions without testing data. The effectiveness of the model was verified using data of an independent specimen and reported creep data on another type of PBM. Furthermore, the performance of the proposed model was compared with an existing full-stage creep model. The performances of the two models in terms of RMSE were quantified and compared. Based on the current results, the following conclusions were drawn.

•
The developed unified phenomenological creep model can describe the full primarysecondary-tertiary creep process under more general conditions of temperature and stress. The effectiveness of the model was validated using both independent and third-party testing data. • The Larson-Miller parameter can be used for predicting the rupture time of PBMs. Combined with the proposed strain rate model, it can be used to predict the creep behavior under more general conditions without testing data on rupture life. • The developed model was compared with an existing reference model. Results show that the developed model is more accurate in terms of root mean squared error. For the testing data used in this study, the proposed model reduces the overall errors by 28.5%. In addition, the proposed model is more reliable for the tertiary creep prediction due to the monotonic strain rate equation.
It is worth mentioning that the proposed model is phenomenological in nature and it cannot explain the detailed creep mechanisms of PBMs. However, it provides a viable means for creep strain and rupture time prediction using one unified model under more general conditions of temperature and stress. At least three sets of testing data on the entire creep process are required to identify the required parameters in the proposed model. Informed Consent Statement: This study did not involve humans.
Data Availability Statement: All data are included within the text.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: PBMs Polymer-bonded composites materials CT Computed tomography RMSE Root mean square error