A General Temperature-Dependent Stress–Strain Constitutive Model for Polymer-Bonded Composite Materials

This study develops a general temperature-dependent stress–strain constitutive model for polymer-bonded composite materials, allowing for the prediction of deformation behaviors under tension and compression in the testing temperature range. Laboratory testing of the material specimens in uniaxial tension and compression at multiple temperatures ranging from −40 ∘C to 75 ∘C is performed. The testing data reveal that the stress–strain response can be divided into two general regimes, namely, a short elastic part followed by the plastic part; therefore, the Ramberg–Osgood relationship is proposed to build the stress–strain constitutive model at a single temperature. By correlating the model parameters with the corresponding temperature using a response surface, a general temperature-dependent stress–strain constitutive model is established. The effectiveness and accuracy of the proposed model are validated using several independent sets of testing data and third-party data. The performance of the proposed model is compared with an existing reference model. The validation and comparison results show that the proposed model has a lower number of parameters and yields smaller relative errors. The proposed constitutive model is further implemented as a user material routine in a finite element package. A simple structural example using the developed user material is presented and its accuracy is verified.


Introduction
Polymer-bonded composites materials (PBMs) have drawn increasing attention in recent years due to their advantage in complex shaping, bi-injection molding, and ease of machining. By compounding different filler materials with polymer, new functional materials can be made, providing a viable technology for on-demand material design and engineering [1]. Applications of PBMs include, but are not limited to, polymerbonded magnetics [2,3], energetics [4,5], aluminum [6], fiber [7,8], and many others [9][10][11]. To utilize those functional materials in critical applications, understanding their mechanical properties and deformation behaviors under different loading and temperature conditions is of great importance in meeting the basic design criteria.
The mechanical properties such as compliance, modulus, toughness, and strength of polymer materials depend largely on strain (rate) and temperature [12,13]. The stress-strain behavior of PBMs is usually that of viscoelastic plastic, which can be approximated as an elastic (linear) regime followed by a plastic (nonlinear) regime [14] in phenomenological view. The deformation is homogeneous in the elastic regime and can be related to the applying stress by Hooke's law [15]. Once the stress of the material reaches a certain stress level, the material will produce irreversible strain and plastic deformation [16].
The maximum stress below which the resulting strain remains proportional to the applying stress is the elastic limit of the polymer materials. However, the filler materials can drastically alter the elastic limit and change the stress-strain behavior in both elastic and plastic regimes [17][18][19]. In addition, the plastic regime of the PBM can also be vastly different to that of the regular polymer materials due to the use of different filler materials.
The need to incorporate temperature effects into the deformation behavior of PBMs has been identified in several experimental studies on a variety of materials, for example, EDC37 in the temperature range of −65 • C to 60 • C by Williamson et al. [20], Rowanex between −60 • C and 60 • C by Walley et al. [21], HTPB at low temperatures by Chen et al. [22], short carbon-fiber-reinforced polyether-ether-ketone (SCFR/PEEK) composites between 20 • C and 235 • C by Zheng et al. [23], epoxy mortar at temperature between 80 • C to 210 • C by Vecchio et al. [24], nanoparticle/epoxy nanocomposites between 23 • C to 53 • C by Unger et al. [25], and many others [26,27]. Several models have been reported to describe the mechanical responses of PBM under different temperatures, including the elastic-plastic fracture-mechanism-based model using the glass-transition temperature (T g ) as a piecewise knot point with a total number of 39 parameters [28], temperature-dependent mesoscale interface-based model [29], and the temperature-dependent visco-hyperelastic model [30]. Models based on the concept of temperature correction factor and its variants have been seen, for example, in [31,32]. Experimental evidence indicates that the mechanical behaviors of PBMs under tension and compression are quite different, and the stress-strain behavior under tension and compression may be described separately [20,21,28]. The aforementioned models in [29][30][31] were established for a single mode, and their applicability to the other mode is unknown. Other models such as the one reported in [28] and its variants requires a dozen parameters, demanding much more testing data for parameter calibration, which may not be realistic in all cases.
The goal of this study is to develop a general temperature-dependent stress-strain constitutive model for PBMs suitable for both tension and compression loading modes with a minimal set of parameters. Compression and tension testing in the temperature range from −40 • C to 75 • C is performed to obtain stress-strain response data. Based on the observation of the stress-strain curves at individual temperatures, a nonlinear model based on the Ramberg-Osgood (RO) relationship is proposed to describe the stress-strain under a given temperature. By correlating the model parameter with temperature, a general temperature-dependent constitutive model can be formulated. The reminder of the study is organized as follows. First, the experimental work is presented in detail and the data in compression and tension are acquired. Next, the constitutive model is developed based on the testing data. Following that, the model validation using both an independent set of testing data and third-party data is performed. The effectiveness of the model is demonstrated, and the performance of the model is compared with a reference model having a comparable number of parameters. The developed constitutive model is further implemented as a user material in a finite element package for structural applications. Finally, conclusions are drawn based on the current results. Figure 1 presents the overall process of experimental testing, model development, validation, and comparison. In the experiment, test specimens were prepared according to the national standard [33] and were arbitrarily divided into two groups, e.g., 80% for model development and 20% for validation. The testing is performed at various temperatures in the range of −40 • C to 75 • C under tension and compression. In the model development part, a model based on the RO relationship is proposed to describe the stress-strain data under each of the temperatures. The temperature dependence on model parameters is established using the response surface model and the final temperature-dependent constitutive model is formulated. In model validation, the validation data (20%) and third party data are used to verify the effectiveness of the proposed model. In model comparison,

Specimens Preparation and Test Conditions
The composition of the PBM used in this study is barium sulfate of 94 wt.% embedded in 6 wt.% fluororubber binder. The conversion volume percentages are 87.6% and 12.4%, respectively. The size of the filler particles is in the range of 0.5 mm and 3 mm. The test specimens were obtained from a vendor and the actual image of the PBM sample is shown in Figure 2a. The PBM was inspected using cone-beam computed tomography (CT), as shown in Figure 2b, to ensure the integrity before testing. In the shown CT image, taken from a piece of the material sample, the light-colored particles are filler materials and the dark-colored portion are the binder materials. After CT inspection, the cylinder specimens for compression (Φ20 mm × 20 mm) and dog-bone specimens for tension (Φ15 mm × 50 mm) were fabricated according to the standard [33], with the dimensions shown in Figure 3a,b, respectively.  A uniform scheme of design of experiments was used to determine the design points of the temperature. In this study, 20 test conditions are designed in the temperature range of −40 • C to 75 • C and the temperature interval is less than 10 • C, totaling a number of 191 test specimens. The temperatures and the number of specimens at the given temperature are summarized in Table 1. Moreover, one validation specimen is selected for every 4 points, excluding the boundary temperature points. Table 1. Test temperature, the number of repeated specimens, and the usage of data. The symbols 'M' and 'V' denote modeling and validation, respectively. Striking lines denote no testing performed at that temperature.

Tension Compression
It is observed from Figure 4 that the ultimate strength decreases significantly as the temperature increases, e.g., from 11.55 MPa at −40 • C to 2.69 MPa at 75 • C. The compression results shown in Figure 5 can be divided into two regimes. As the strain increases, the stress increases monotonically until the ultimate strength is reached, following by the postultimate stage before fracture. With the further increase of the strain, the stress decreases gradually. For compression results, the ultimate strength decreases from 42.41 MPa at −40 • C to 12.64 MPa at 75 • C. The maximum resistance of a material to external forces is called the stress at break [34]. The ultimate strength in tension and compression indicates the failure. Therefore, the data ranging from the initial state to the stress at break are used for the model development. The main transition temperature for the used binder material is T g = 80 • C, which is consistent with that of a similar material [35]. The stress-strain behavior at a higher temperature that is close to T g may be influenced by this transition temperature. This can be observed from Figures 4 and 5, where at higher temperatures, the stress-strain curves start to deviate much more than those at lower temperatures. As the proposed constitutive model is empirical in nature, the developed model does not deal with the mechanism related to the transition temperature.

Constitutive Model Development
The stress-strain model under single temperature is proposed first, and the temperature dependence of the model parameters is established using the tension and compression data. After that, the general temperature-dependent stress-strain constitutive model is obtained.

Stress-Strain Model under Single Temperature
The stress-strain curves presented in Figures 4 and 5 indicate that the total response consists of a short linear part and a nonlinear part. Therefore, the total strain ε can be written as the linear accumulation of the two parts such as where ε e is linear (elastic) strain and ε p is nonlinear (plastic) strain. For the linear strain, it can be described by the stress σ and Young's modulus E as (σ/E). For the nonlinear part, a general exponential form as in the Ramberg-Osgood relationship [36] can be used. Substituting the linear and nonlinear parts, for monotonic loading, to determine the adequacy of Equation (2) in describing the testing data at a single temperature. The testing data of T = −40 • C are arbitrarily chosen to fit the model parameters, and the fitted model parameters are presented in Table 2 for both tension and compression cases. Using the fitting parameters, the results of the mean fit and the testing data are presented in Figure 6 for comparison purposes. It can be observed that the model can capture the linear and nonlinear portions of the testing data fairly well in the given strain range for both tension and compression cases, indicating the effectiveness of Equation (2) for correlating the stress-strain behavior of the PBM under investigation.

Model Parameters under Different Temperatures
The nonlinear least squares estimator is used to obtain the parameters (E, K, n) in Equation (2) for tension testing data of each of the specimens at each of the temperatures. For convenience, the parameter ln K instead of K is used. The resulting parameters vs. temperatures are presented in the Figure 7. Every solid dot in Figure 7 is the fitted parameter associated with one specimen. From the results of Figure 7, it can be found that the variation of each of the three parameters with temperatures can loosely be divided into two linear regimes. For parameters ln K and n, the change points of the two regimes are about −20 • C, indicated by the vertical dash lines in Figure 7a,b. The two parameters decrease linearly from −40 • C to the change point, and increase linearly from that point to the end. For parameter E, the change point is also −20 • C but the variation in the two regimes is reversed compared to that of the other two parameters. The scattering in the fitted parameters shown in Figure 7 is caused by the inherent uncertainty in material properties and stochastic nature of the filler particle distribution. -50

Temperature Dependence Modeling
Based on the above discussion, it can be observed that the temperature dependence can be described using the piecewise linear relationship; therefore, a piecewise linear model is used to correlate the variation of the three parameters with temperatures.
where α i , β i , and γ i (i = 1, . . . , 3) are fitting coefficients without unit; T is the temperature in degrees centigrade; and T ki (i = 1, . . . , 3) are the knot points of the piecewise linear function. The term H(x) is the Heaviside step function defined as According to Figure 7, T ki (i = 1, . . . , 3) are determined as −20 • C. Using the parameter-fitting results presented in Figure 7 and Equation (3), the temperature dependence fitting coefficients α, β, and γ are obtained using the ordinary least square estimator as (5) Figure 8 presents the temperature dependence fitting results, where the mean fit (solid lines) is shown. It can be seen that the mean fit of the piecewise linear model in Equation (3) can roughly describe the temperature dependence of the model parameters ln K, n, and E.  (5), the final temperature-dependent stress-strain constitutive model under tension is obtained as

Substituting Equation (3) into Equation (2) with parameters in Equation
where the temperature terms are Tension : The temperature-dependent stress-strain model (Equation (6)) is used to predict the strain given the temperature and stress. The overall residuals on all the tension data are evaluated and presented in Figure 9. The value corresponding to the sum of squared errors (SSE) is 3.0686 × 10 −4 using Equation (8). The SSE is given as where y i is the actual value,ŷ i is the model prediction, and i = 1, . . . , N represents the index of a total number of N data points. The standard deviation of the residuals is estimated as 5.133 × 10 −5 .

Model Parameters under Different Temperatures
The same fitting and processing method in Section 3.2.1 is applied to the compression data, and the resulting parameters in Equation (2) under different temperatures are presented in Figure 10. Similar to the tension results, the variation of parameters with the temperatures under compression are also piecewise linear. Both ln K and n are linear proportional to the temperature in the range of −40-40 • C. In the range of 40-75 • C, the two parameters are inversely proportional to the temperature. The parameter E shows the same behavior but with a different change point of −5 • C. Therefore, for compression data, the same equation, Equation (3), is used with T k1 = T k2 = 40 • C and T k3 = −5 • C.

Temperature Dependence Modeling
Using the parameter-fitting results presented in Figure 10 and Equation (3), the temperature dependence fitting coefficients α, β, and γ are obtained using the ordinary least squares estimator as Substituting Equation (3) into Equation (2) with parameters in Equation (9), the final temperature-dependent stress-strain constitutive model under compression is obtained as where the temperature terms are  Figure 11 presents the temperature dependence fitting results, where the mean fits (solid lines) are shown. It can be observed that the temperature dependence can be described using the piecewise linear functions. Figure 12 presents the residuals of the model fitting for compression data. The SSE in this case is 0.07323 and the standard deviation of the residuals is 0.001089.

Model Validations and Comparisons
Independent sets of testing data are used to validate the performance of the proposed model. Third-party data published for a PBM are further used to validate the generality of the model. To compare the proposed model with existing models, a model with a comparable number of parameters is adopted and results from the two models are compared. In the validation and comparison, the following two general measures, namely, the mean absolute error (MAE) and the root mean square error (RMSE), are used.
where y i ,ŷ i , i, and N are defined as in Equation (8).

Validation Using Independent Sets of Testing Data
The data from the validation specimens (labeled as 'V' in Table 1) are not used for model development and are used for validation. It is noted that both tension validation specimens and compression validation specimens are available. We utilize data associated with the validation specimens. Prediction results of the developed model and the actual testing values are compared. In addition, 2 groups of modeling data near the piecewise point are added to compare with the model fitting results. Figure 13 presents the model prediction results for tension specimens compared with testing results. Figure 14 shows the model prediction results for compression specimens compared with testing results. In both modes, a close agreement between the prediction and testing data can be observed.
The maximum MAE and RMSE are 6.200 × 10 −5 and 8.625 × 10 −5 at 25 • C, respectively, for tension data. For compression data, the maximum MAE and RMSE are 0.001741 and 0.002460 at 10 • C, respectively. The prediction performance for each of the temperatures are detailed in Table 3.

Validation Using Third-Party Data
To further verify the generality of the proposed model, third-party testing data on PBX 9502 from [28] consisted of 5 sets of tension data and 5 sets of compression data (at temperatures −52 • C, −20 • C, 20 • C, 50 • C, and 72 • C) are used to validate the generality of the proposed model.
It should be noted that the results of model parameters in Equations (7) and (11) cannot be directly used for prediction purposes as those parameters are material specific. The reported data on PBX 9502 are used to fit the developed model using Equation (2) for tension and compression data at each of the temperatures, and the temperature dependence for both modes are obtained using Equation (3). The same least squares fitting scheme is used. The prediction results of the proposed model and the actual testing data at the five temperatures are compared and shown in Figure 15. Under tension, the maximum MAE and RMSE are 5.4472 × 10 −5 and 7.8292 × 10 −5 at 50 • C, respectively. For compression, the maximum MAE is 0.001040 at 74 • C and the maximum RMSE is 0.001184 at 20 • C. For both cases, the model can capture the behaviors of the stress-strain variations fairly well.

Model Comparisons
To further investigate the performance of the proposed model, a stress-strain constitutive model reported by Chang et al. [37] is adopted for comparison purposes. This model is formed by appending a quadratic term and a cubic term of the strain to the regular linear model, as in Equation (14): where the parameters E 0 , C 1 , and C 2 are fitting parameters. It is worth mentioning that other reported models can also be adopted for comparisons; however, this model is chosen due to the fact that the model has a comparative number of parameters and it requires only stress-strain data for fitting. Other models reviewed in the introduction either contain too many parameters or require special testing data.

Tension
The same amount of tension data is used to fit model parameters using Equation (14) for each of the temperatures. The fitted model parameters (E 0 , C 1 , C 2 ) are further correlated with temperatures to obtain the following temperature dependence equations in Equation (15): The two models with established model parameters are used to predict the validation data (those labeled as 'V' in Table 1), and the results are shown in Figure 16. In general, both models can describe the data satisfactorily; however, for data at 50 • C, the proposed model outperforms the reference model. The errors in terms of MAE and RMSE of the two models are presented in Table 4. From Table 4, it can be seen that the proposed model yields smaller errors for all data sets except the one associated with 0 • C.

Compression
The same amount of data under compression are used to fit the reference model, resulting in the following temperature dependence functions in Equation (4).
The comparisons for compression data between the two models are made following the same procedure described above and are shown in Figure 17, with the error results compared in Table 5. In general, the proposed model can yield improved prediction results, particularly for the tail regions of the stress-strain data. The detailed error results in Table 5 indicate that the proposed model outperforms the reference model at all temperatures except 10 • C. It should be noted that other probabilistic methods such as Bayesian model assessment [38][39][40] can also be employed to compare the performance of the models.

A User Material Implementation of the Model
To use the proposed temperature-dependent stress-strain model for structural applications, a user material subroutine for the proposed model is developed for finite element analysis. The dimensional of the structural model is the same as the actual specimens. The structural models for the tension and compression specimens are created and meshed using quadratic hexahedron elements with an average size of 2 mm.
For the tension specimen, a constant stress of 2 MPa is applied to the top face with the bottom face fixed in all degrees of freedom, representing a tension test condition. For the compression specimen, a constant stress of −20 MPa is applied on the top face, and the same boundary conditions are used for the bottom face. Figure 18 demonstrates the resulting normal stresses after applying the load at −40 • C, which are shown for the tension and compression specimens in Figure 18a,b, respectively. Under tension, it can be found that the stress on the top face is about 2 MPa, which is the same as the applied load, and the maximum stress occurs in the middle section of the structural part. Under compression, the maximum stress appears at the fixed position. From the top face, the stress gradually decreases along the UX direction. The stress on the top face is the same as the applied load of 20 MPa. For both cases, the resulting stresses are consistent with the actually applied loads.  Figure 19 presents the FEA results using the specimen structural models in Figure 18 with the developed user material routine for temperatures of −40 • C, 20 • C, and 75 • C and compares them with the theoretical results of the proposed temperature-dependent stress-strain model. Both tension and compression results are compared. The maximum relative error (RE) is 0.01239, occurring at −40 • C for the tension case, and the maximum RE is 0.01537 at 75 • C for the compression case, indicating the effectiveness of the developed user material routine for structural applications.

Conclusions
A general temperature-dependent stress-strain constitutive model for polymer-bonded composite materials is proposed, developed, validated, and compared. Experimental testing data were acquired in the temperature range of −40-75 • C for both tension and compression. The model under a given temperature is proposed based the linear accumulation of an elastic term and a plastic term, resembling the Ramberg-Osgood relationship. The model parameters are obtained by fitting the testing data at each of the temperatures, and the temperature dependence of the model parameters can be established using the response surface method. The proposed model is validated using independent sets of testing data and third-party data, and is further compared with an existing model with a comparable number of parameters. A user material routine of the proposed model is implemented and verified for structural applications. Based on the current results, the following conclusions are drawn.

•
The proposed method can unify the temperature effect by modeling the temperature dependence of the constitutive model parameters, thus providing an alternative to the existing temperature correction factor-based methods. Furthermore, such a treatment allows a minimal number of fitting parameters compared to the existing models. • The proposed method provides a general approach to model the stress-strain behaviors of PBMs. The basic assumption that the stress-strain response is an accumulation of an elastic term and a power-law plastic term is observed in many other materials such as various metals. The effectiveness and generality of the model are validated and demonstrated by two independent preponderate data sets in the temperature range −40-75 • C.
It is worth mentioning that the developed model is a semiempirical data-driven model, and it does not deal with the microscopic mechanism. The dynamic mechanical analysis on the material can be carried out to enhance the understanding of the mechanism, and SEM imaging may be used to reveal the interface interaction between the filler particles and the binder materials. The constitutive model established in this study does not incorporate the strain rate effect. The strain-rate-induced hardening phenomenon has been reported in several studies, e.g., [41,42]. The incorporation of the strain rate effect to the proposed model will be studied in future work. The proposed model considers the PBM as a homogeneous material in macroscopic perspective. The local variations of the material properties due to the microscopic effects of filler materials and polymer binder interface [43][44][45] may be explored in the future. In addition, the effectiveness of the model is verified using the presented experimental data acquired in this study and a third-party published data. Due to the limit on currently accessible data, the verification of the developed model using other PBMs will be investigated in the future. The current modeling adopts a deterministic approach, and a probabilistic approach can also be employed for reliability analysis [46][47][48]. 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: