Comparison and review of classical and machine learning-based constitutive models for polymers used in aeronautical thermoplastic composites

2.1.1 Rubber-like models

This kind of models are developed through the elastic theory of rubber, so they could be called rubber-like models. Studies of rubber elasticity began as early as the early twentieth century, and James and Guth [20] discussed an idea of effective internal pressure, and the typical stress–strain curve for rubber is displayed in Figure 3. They assumed that the polymer chains move through each other only through crosslinks, and the whole system is prevented from collapsing by the assumption of repulsive forces which can generate a bulk modulus. In addition, a simplified model for the bulk rubber was proposed, which consists of a network of the idealized flexible chains extending through the material and fluid filling it. Moreover, the bounding surfaces are in equilibrium under all the forces acting on them (internal pressure, pull of the molecular network, and any external forces). Because of the changing entropy, a good agreement was obtained by comparing the stress–strain curves for bulk rubber at a constant temperature. The general theory proposed by James and Guth [20] provides a basis for the treatment of other physical properties of stretched rubber.

Figure 3

Typical stress–strain curve for rubber [20].

Ball et al. [21] indicated that the assumption by James and Guth [20] ignored the repulsive forces which generate a bulk modulus. They proposed that entanglements can be simulated by links that make a sliding contact between polymer networks. The contribution of entanglement to the free energy of shear is given in Eq. (1), where λ _i is the Cartesian extension ratio and η is a measure of the freedom of a link to slide compared with the freedom of movement of a chain.

Edwards and Vilgis [18] proposed an elasticity theory of rubber based on the concept of entanglements in 1986, which is an improvement of the work presented by Ball et al. [21]. They found that the molecular mechanism of stretching is dominated by the slippage of chains for small deformation, and the hardening of the rubber at high deformation attributes to the inextensibility as described by the tube concept. The constitutive model is able to explain the deformation behavior of rubbers over the total range of deformation, by comparing the free energy of deformation with the experiment. The free energy is shown in Eq. (2), where α is a measure of the inextensibility and η of the slippage, N _c is the number of crosslinks, and N _s denotes the number of slip links.

The Edwards–Vilgis model was the most representative one in this type of approach, and it was also known as the network model. It was initially applied to rubbers, then extended and applied to amorphous glassy polymers or semi-crystalline polymers by other researchers. Van Ruiten et al. [22] applied the Edwards–Vilgis model for the analysis of melt-spun PA 4.6 fibers in 2001. They tested the maximum attainable tenacity of drawn yarns under given drawing conditions, and compared the experimental results with the predictions of the Edwards–Vilgis model. It was shown that the Edwards–Vilgis rubber-elastic model can describe the network deformational behavior of the as-spun yarns over a wide range of draw ratios. Based on the rubber elastic model, Sweeney et al. [23] suggested a large deformation and rate-dependent model in 2002. In their theory, the rigid spheres are embedded to introduce strain concentration, which is similar to that caused by hard crystalline regions. Besides, the dependence on time and rate was introduced via the shear stress-driven diminution of the sphere radii. Then, the model was applied to the high temperature stretching of PE through ABAQUS software (a finite element (FE) solver). Constitutive models introduced in this section are listed in Table A1 of the Appendix.

2.1.2 Models based on Haward–Thackray theory

This type of model takes the theory of Ree and Eyring [24] as the foundation. Ree and Eyring [24] presented a general relationship between the yield stress, the strain rate, and the temperature applicable to polymers as early as 1955, and Ree and Eyring’s theory can account for the relaxation process of viscous flow. Since then, researchers have come up with models based on this theory.

Haward and Thackray [17] represented the interactions between molecules through springs and dashpots in 1968. This model is intended to provide a semi-empirical approach to the large plastic deformations at the yield point in low-temperature where the temperature is below T _g. Nevertheless, it is not suitable for describing the creep behavior at extensions of up to 5% or nonlinearity in the initial Hookean modulus. Haward–Thackray model contains a Hookean spring, an Eyring dashpot, and a Langevin spring with an ultimate limiting network strain as shown in Figure 4. The above Hookean spring represents the constant Hookean modulus, which describes the linear elastic deformation of the amorphous phase. The Eyring dashpot represents the VP character, which describes the rate-dependent macroscopic yield deformation of the amorphous phase. The Langevin spring represents the limited elastic extensibility, which describes the strain hardening result from the change in configuration entropy caused by the orientation rearrangement of molecular chain entanglement. The limited elastic extensibility is based either on the first-order process model, or the conventional Langevin formula as used to describe the highly elastic extension of rubbers. Practice shows that the model can be used to characterize the intermolecular elastic interaction, strain-rate dependent yield behavior, and strain-strengthening behavior in the post-yield stage of amorphous polymers.

Figure 4

Haward–Thackray model system [17].

Haward–Thackray theory is a milestone and many subsequent models are inspired by it to varying degrees. On this basis, a large number of constitutive models are proposed. We introduce the development profile with three typical types of models and their derivatives: the BPA model, the glass-rubber (GR) model, and the Eindhoven glassy polymer (EGP) model.

2.1.2.1 BPA model and its derivatives

Among numerous research groups, Boyce et al. [25] have made an outstanding contribution. Boyce et al. [25] proposed a three-dimensional constitutive model for glassy polymer on the basis of the Haward–Thackray model in 1988, also known as the BPA model. Linear spring, VP dashpot, and Langevin spring are contained in the BPA model as shown in Figure 5, and the linear spring and VP dashpot are connected in series to construct the Maxwell element. BPA model assumes that the deformation resistance of amorphous polymers may be decomposed into intermolecular deformation resistance and entropy deformation resistance. The intermolecular deformation resistance is related to the rotation of molecular chain segments and determines the elastic deformation before yield. The entropy deformation resistance is related to the rearrangement of molecular chain segments, and will increase due to the continuous rearrangement of molecular chain segments during the post-yield process. Strain hardening is controlled by entropy deformation resistance, while strain softening is controlled by both intermolecular deformation resistance and entropy deformation resistance. BPA model can achieve good results by comparing with the experimental data of PMMA.

Figure 5

BPA model [25].

2.2.1 Elasto-plastic models

G’sell et al. [34,35] proposed a phenomenological model through a series of tests of polyvinyl chloride (PVC) and high-density polyethylene (HDPE) in 1979. The stress can be expressed as a function of strain, strain rate, and temperature in Eq. (3).

where σ, ε, ε ̇ , and T denote stress, strain, strain rate, and temperature, respectively, and K, W, and m are the model parameters. The G’sell–Jonas model is suitable for the deformation prediction of materials at a low strain rate (10⁻¹–10⁻⁴ s⁻¹). Johnson and Cook [36,37] proposed a famous Johnson–Cook model (JC model) when studying the deformation of metal under high strain rate and high temperature in 1985 as shown in Eq. (4).

where σ, ε p , ε ̇ p , ε ̇ 0 , T, T _m, and T _ref signify stress, plastic strain, plastic strain rate, reference strain rate, temperature, melting temperature, and reference temperature, respectively, and there are only five model parameters in the JC model (A, B, C, n, and m). JC model is a typical model, and it reflects the behavior of elasticity, strain hardening, and strain rate hardening for metals. However, researchers subsequently found that the JC model was also suitable for predicting the deformation of semi-crystalline polymers. Garcia-Gonzalez et al. [38] applied the JC model to the simulation of impact for PEEK and Ti6Al4V titanium alloy. In addition, the JC model was constantly extended in order to get better performance. For example, Chen et al. [39,40] modified the original JC model into the following form shown in Eq. (5), so that the modified one has better capability to predict the flow behavior at elevated temperature conditions.

Neither of G’sell–Jonas model and JC model mentioned above can accurately describe the strain softening deformation of the polymers. Thus, Duan et al. [41] proposed a homogeneity model suitable for glassy and semi-crystalline polymers based on G’sell–Jonas model, JC model, Matsuoka model, and Brook model, named DSGZ model. DSGZ model can describe the characteristics of elasticity, yielding, strain hardening, and strain softening by taking strain rate and temperature into account, as shown in Eq. (6).

where C ₁, C ₂, C ₃, C ₄, K, a, m, and α are model parameters, and they are easy to be calculated according to the strain–stress relationship. Then, Duan et al. [42,43] applied the DSGZ model to predict the deformation under the impact load of ABS, PC, and polybutylene-terephthalates (PBT). Further, the DSGZ model has been continuously improved and applied to ABS, PBT, PC, PMMA, PA, and PEEK by researchers [12,44,45,46,47].

In order to further evaluate the performance of phenomenological models, the deformation of PMMA and PEEK at a low strain rate was simulated by applying the DSGZ model by the authors. The experimental data of PMMA was taken from the article of Duan et al. [41], and the experimental data of PEEK was taken from the research of Chang et al. [48]. Particle swarm optimization was applied to fit the parameters of the DSGZ model according to the experimental data. By comparing the results of prediction with the experiment (as shown in Figure 6), it was found that the DSGZ model has the ability to capture the behavior of polymers at different temperatures. Moreover, the nonlinear elastic deformation, yielding, strain softening, and strain hardening at low strain rates were captured well with a reasonable fitting of parameters. Although the prediction of the DSGZ model may not be guaranteed at a high strain rate, the phenomenological constitutive models can always find appropriate forms to obtain a satisfactory prediction.

Figure 6

Comparison of true stress–strain curves between the prediction of DSGZ model with experimental results, left is for PMMA and right is for PEEK.

Mulliken and Boyce [49] proposed an M–B model on the basis of the BPA model in 2006. They decomposed the intermolecular resistance into two rate-dependent resistances as shown in Figure 7. In other words, network A in the BPA model was decomposed into two parallel Maxwell elements which are related to glass transition (α) and secondary molecular motions (β), respectively. In this way, the behavior of polymers in low strain rate can be reflected by α, and the behavior of polymers in high strain rate is determined by both α and β. The M–B model can not only capture the transition in the yield behavior, but also accurately predict the post-yield and large strain behavior over a wide range of temperatures and strain rates. The total stresses in the polymer are given as the tensorial sum of the intermolecular stresses of α and β and the network (back) stress. Moreover, expressions of three stress components are presented in Eq. (7). It has been indicated that the M–B model can accurately predict the behavior of materials in the strain rate range of 10⁻⁴–10³ s⁻¹.

Figure 7

M–B model [49].

The M–B model is capable of predicting the deformation trend of polymers at low, moderate, and high strain rates. Figure 8 shows that it is convenient to use M–B model for analysis of many polymers such as PC (shown in Figure 8 (a)) and PMMA (shown in Figure 8 (b)), because the deformation response is divided into α and β. However, it was found that the M–B model ignored the viscosity effect before yield. In addition, the strain softening stage predicted by the M–B model usually deviate from the experimental results due to the early arrival of yield. The M–B model can only predict the trend of strain hardening at a low strain rate, but the prediction becomes inaccurate with the increase in the strain rate. In addition, the M–B model is an isothermal model, and it cannot capture the post-yield thermal softening. Other researchers have drawn similar conclusions [50].

Figure 8

Comparison of true stress-strain curves between the predictions of M–B model and experimental results, where (a) shows the contrast for PC and (b) shows the contrast for PMMA [49].

By combining the M–B model, G’sell model, and DSGZ model, Wang et al. [50,51] proposed an adiabatic phenomenological constitutive model to predict the mechanical behavior of PC at various strain rates and temperatures. Deformations of the low strain rate and the high strain rate were reflected by the α-component and the β-component, respectively. In addition, they used the model to simulate the split Hopkinson pressure bar testing and the falling weight impact testing by explicit user material subroutine (VUMAT). Moreover, results of the simulations agree well with the experimental data. In addition, some phenomenological models have been also developed for 3D-printed polymers [52].

Varghese and Batra [53] modified the M–B model by introducing temperature and strain rate-dependent elastic moduli and two internal variables at the high strain rate. Further, Safari et al. [54,55] indicated that α and β alone in network A were insufficient under the high strain rate, so γ transition as shown in Figure 9 was added to reflect the behavior of materials in the higher strain rate. Moreover, their studies have shown that the modified model can predict the thermomechanical behavior of polymer when the strain rate is over 10,000 s⁻¹. However, Safari et al. [54,55] indicated that due to adiabatic condition of high strain rate deformations, the modified model requires higher strain rates to account for temperature changes.

Figure 9

Safari model [54].

2.2.2 VE models

Both reversible and irreversible deformation is usually included for polymers, and some researchers argue that the reversible deformation of polymers should be concerned. Nonlinear VE characteristics will thus appear in such types of constitutive models. Schapery [56] proposed that the behavior of some steels and polymers such as creep and relaxation can be reflected by linear equations, nonlinear equations, and some other specific equations which are summarized from a good deal of experiments. In fact, Schapery [56] attempted to build some connections between experimental phenomena and some constitutive equations, based on a framework of irreversible thermodynamics that are similar to the Boltzmann superposition integral form of linear theory. They suggested the core idea of phenomenological constitutive models and lay a foundation for the subsequent development of phenomenological constitutive models. Khan et al. [57] proposed a one-dimensional (1D) phenomenological constitutive model with a semi-empirical modification based on experimental observations, to capture the complex and highly nonlinear finite thermo-mechanical behaviors of VE polymers. Moreover, the model was built based on the infinitesimal linear theory. It was demonstrated that the model was accurate in predicting the deformation of polymers over a wide range of strain rates at room temperature or near the T _g. Chang et al. [48] proposed a phenomenological nonlinear VE constitutive model to characterize the stress–strain of PEEK before yielding with considering temperature and strain rate dependence. The VE model consists of a cubic nonlinear spring and a linear Maxwell element connected in parallel, and the form of this model as shown in Eq. (8) is similar to the subsequent Z-W-T model [58,59].

Wang et al. [58,59] proposed a Z-W-T model to represent the nonlinear VE deformation of polymers in 1991, and then they applied this model for the impact investigation in bird strikes on windshields of the high-speed aircraft. The expression is shown in Eq. (9).

where E ₁, θ ₁, and E ₂, θ ₂ signify the elastic modulus and relaxation time of VE response at low or high strain rate, and E ₀, α, and β are nonlinear elastic modulus except for the VE response. In the early stage, the Z-W-T model was only suitable for polymers at room temperature, because the temperature was not considered in the model. Therefore, Wang et al. [60] extended the Z-W-T model by introducing the effect of temperature, and simulated the dynamic response of PMMA windshield against bird strike accurately. Dar et al. [61] also implemented the Z-W-T model to describe the mechanical behavior of PMMA. Nevertheless, the Z-W-T model is adaptable to simulate the VE response but not the VP response after yielding.

2.2.3 VP models

Some researchers pay attention to the irreversible deformation, and suggested a series of VP constitutive models. A series of models were improved based on over-stress model (VBO). Colak [62] proposed a modified viscoplasticity theory based on the over-stress model, which contains two tensor valued state variables, the equilibrium, and kinematic stresses and two scalars valued state variables, drag, and isotropic stresses. They compared the numerical results with the experimental data, and the tests of polyphenylene oxide were performed at different stresses above and below the yield point. It is shown that nonlinear rate sensitivity, nonlinear unloading, creep, and recovery at zero stress were simulated by the proposed model. Besides, Ghorbel [63], Drozdov [64], Khan and Yeakle [65], Dusunceli and Colak [66] developed the corresponding models based on the VBO model for analysis of the polymers such as PA, PC, PET, and PP. Bardenhagen et al. [67] presented a general framework to develop the constitutive models of polymeric materials in the VP regime, and a VP constitutive model of three-dimensional finite deformation. Strain-rate dependence, stress relaxation, and creep phenomena can be reflected in the proposed model. Drozdov and Christiansen [68] presented a VP constitutive model for isothermal three-dimensional cyclic deformations with small strains of semi-crystalline polymers. Moreover, multiple inelastic deformations appear in the model with the 1D spring and dashpot construction. Although there are 15 material parameters in the proposed model, they can be determined step by step by matching appropriate intervals of a stress–strain curve at loading and retraction. Good prediction can be obtained, when the constitutive model was applied to predict the VP behavior of HDPE with a strain less than 0.1.

2.2.4 Hyperelastic models

Neo-Hooken model [69] is a classic statistical thermodynamic hyperelastic model, and it is suitable for the deformation of rubber-like materials. The Neo-Hooken model is expressed for its Helmholtz free energy per unit reference volume as shown in Eq. (10), which is independent of the temperature.

where Ψ ( I 1 , J ) is the Helmholtz free energy, I ₁ is the first invariant of strain tensor, and J is the determinant of deformation gradient tensor. In addition, Mooney–Rivlin model [70] is also a classic hyperelastic constitutive model.

Wang and Guth [71] established an orthogonal non-Gaussian three-chain network constitutive model to describe the large deformation of polymer chain segments by statistical methods in 1952. Flory and Rehner [72] proposed a four-chain constitutive model of regular tetrahedra to predict the strain hardening deformation of hyperelastic materials. Arruda et al. [73] proposed a fully three-dimensional eight-chain thermodynamically consistent constitutive model in 1995, and a thermo-mechanically coupled FE analysis for the large deformation of glassy polymer was carried out. The inelastic deformation of strain hardening in this model is not dissipative, and it is to be stored as internal back stress. The strain rate and temperature-dependence of initial yield are included in the model as well as the temperature dependence of evolving anisotropy and its associated strain hardening. Wu and van der Giessen [74] also proposed an eight-chain constitutive model, which features strain-rate, pressure, and temperature-dependent yield, softening immediately after yield, and subsequent orientational hardening with further plastic deformation. Sweeney and Ward [75] conducted the experiments near and above the glass transition temperature (T _g) of PVC in 1995, in which PVC is an amorphous polymer containing a small amount of crystallization. In addition, they introduced strain-rate into Edwards–Vilgis model, and investigated the strain-rate dependence of PVC. Billon [76] proposed a 1D visco-hyperelastic model, and validated it by using a rich and rigorous database of PMMA above T _g. Moreover, they extended the Edwards–Vilgis model to a general time-dependent constitutive model. The inelastic phenomena in this model were reflected by an evolution of internal variables related to the alteration of microstructure, that induces changes in parameters in constitutive model and dissipation of energy. This approach worked well in modeling time-dependent behavior of polymers over a wide range of temperatures and strain rates. Maurel-Pantel et al. [77] extended the model proposed by Billon [76] to three-dimensional thermodynamically consistent constitutive equations, and proposed a visco-hyperelastic constitutive model in 2015. On the basis of Edward–Vilgis theory, they also represented the degree of mobility of entanglement points through the introduction of an evolution equation for the internal state variable. The thermomechanical model was applied to a semi-crystalline polyamide polymer (PA66), and the predicted deformation and temperature were in good agreement with the experimental results under tension and shear conditions. Gehring et al. [78] modified the network theory of Billon [76] and Maurel-Pantel et al. [77] in 2016. For the behavior of amorphous and semi-crystalline PET, they considered several contributions of the microstructure’s rearrangement (disentanglement and loss of connectivity) within the large deformation formalism by the thermodynamic framework. In addition, the model considers microstructure at a mesoscopic level through the description of an equivalent network evolving with internal state variables. Good agreement was obtained by comparing with the results of the nonmonotonic tensile test. Federico et al. [79] used the theory proposed by Gehring et al. [78] to investigate the mechanical behavior of amorphous PMMA with different molecular weights in 2020. The model can account for the elastic contribution of an equivalent network, which experiences inelastic mechanisms coming from the evolution of internal state variables. Therefore, the model [78] possesses the capability to capture the mechanical response of materials at different temperatures and strain rates through the VE and rubbery regimes. In addition, there are still lots of models based on such strategies [19,80–84].

2.2.5 Multi-mode models

However, some features have been ignored logically when constructing the above constitutive models. Different polymers usually exhibit different mechanical properties containing more than one pattern, due to the different degrees of crystallization. The properties of polymers would be the most important features to be considered for the study of constitutive models. Researchers hence attempted to couple the VE character and VP character or more characteristics in their models to obtain more reasonable results.

Anand et al. [85,86] developed an elastic-viscoplastic thermodynamically consistent constitutive model for the amorphous polymers with strain rate and temperature-dependent large-deformation behaviors. Moreover, they implemented their model for simulation of large-strain compression experiments for PMMA, PC, and a cyclo-olefin polymer (Zeonex-690R). It was validated that the model can reflect the strain rate and temperature-dependent yielding, strain-hardening, unloading response, and the temperature transforming. Frank and Brockman [87] presented a constitutive model which combines nonlinear VE character and VP character of polymers for multi-axial isotropic deformation in 2001. Several traditional constitutive models expressed by some simple empirical relationships were included in their model such as linear elasticity, linear viscoelasticity, and plasticity, so that the time-dependent and nonlinear behaviors can be captured. The proposed model was applied to simulate the impact behavior of PC, and the feasibility and accuracy of the model are demonstrated.

Miled et al. [88] proposed a coupled VE–VP constitutive model for homogeneous and isotropic polymers. Moreover, fully implicit integration, a two-step return mapping corrector strategy, and a consistent tangent operator were used in their work. The model is restricted to the regime of small perturbations (small strains, displacements, and rotations), and the VE part of the model is linear and isothermal. The total strain in this model was assumed to be a sum of VE strain and VP strain, and this behavior is a common strategy for multi-mode models. The model was implemented in a three-stage shear test and the simulation of HDPE through ABAQUS software via a user-defined routine (UMAT), and the model was proved to be effective. Yu et al. [89] proposed a nonlinear VE–VP cyclic constitutive model based on the tests of multilevel loading–unloading recovery, creep-recovery, and cyclic tension–compression/tension ones. They extended Schapery’s nonlinear VE model [56] as the VE part of the model. Moreover, as for the accumulation of irrecoverable VP strain produced during cyclic loading, they adopted the Ohno–Abdel-Karim’s nonlinear kinematic hardening rule [90]. This model was verified by comparison of the prediction with the tests of PC. Other researchers [13,91–98] have also made significant contributions in this regard.

2.2.6 Constitutive models with damage

It is out of their depth for the above models to describe the damage or crack during the deformation, even though most of the models have reflected the viscous, elastic, and plastic behavior of polymers well. Therefore, researchers take the damage to materials into account when constructing the constitutive models.

Zairi et al. [99,100] proposed a constitutive model to describe the elasto-viscoplastic damage behaviors of polymers according to both micromechanical (Gurson potential) and phenomenological (modified Bodner–Partom model) models, which contains hydrostatic and failure evolution terms. They investigated the macro-mechanical response and damage micro-mechanism by failure growth in rubber-toughened glassy polymer by the simulation of RTPMMA and HIPS, and a good agreement was obtained. An exploratory approach is provided by this model to capture the damage quantitatively, while some improvements are still essential. Tehrani and Abu Al-Rub [101] proposed a nonlinear VE, VP, and visco-damage constitutive model, and applied the proposed model for the damage evolution simulation of PMMA embedded with silicate nanoclay particles. They introduced a damage evolution variable which is a function of actual stress, hydrostatic stress, total strain, strain rate, temperature, and damage process. A series of simulations were performed, and the model was proved to be reasonable. Balieu et al. [102] proposed a phenomenological non-associated elasto-viscoplastic model coupled with damage in the finite strain framework, which was used to simulate the behavior of a 20% mineral-filled semi-crystalline polymer for a large strain rate range. Perzyna-type VP formulation was used to represent the nonlinear strain-stress response, and the specific non-associated VP potential was used to capture the volume change of the material under tensile and compression loading. In addition, hydrostatic pressure was introduced to reflect the yield stress coupled with the damage model. They implemented the model into FE simulation through UMAT, and the result obtained were in good agreement with the experimental data.

Krairi and Doghri [103] proposed a constitutive model by coupling viscoelasticity, viscoplasticity, and ductile damage within the framework of irreversible thermodynamics. They simulated the time and strain rate dependencies, the Bauschinger effect, and ductile damage through a series of simple expressions. Ductile damage evolution is related to VP strains, and it evolves only with changes in VP strains. The model was validated by tests for different polymers under various loadings. Praud et al. [104] also proposed a VE, VP constitutive model considering ductile damage within the framework of thermodynamics, but the VP behavior of their model is different from that of Krairi and Doghri [103]. In the model proposed by Praud et al. [104], the VE part is described with the help of a series of Kelvin–Voigt branches, while the VE properties of the model of Krairi and Doghri [103] are expressed through Prony series through an integral formulation. Khaleghi et al. [105] suggested a thermodynamically consistent damage model to predict failure of glassy polymers based on the EGP multi-mode model. In their model, the effects of plastic deformation and hydrostatic stress on damage evolution were considered. In addition, Abu Al-Rub et al. [106], Seidel et al. [107], Voyiadjis et al. [108], and Cayzac et al. [109] also made great contributions in this aspect. In addition, some issues had also attracted the attention of investigators, such as the influence of moisture [110], heating [111], finite strains [112], change in crystallization [113], hardening behavior [114 and softening behavior [115]. Phenomenological models can be constructed by combining the physical mechanisms and the results of tests, so phenomenological model can usually predict the deformation well as the model is constructed correctly. In addition, phenomenological constitutive models introduced in this section are listed in Table A3 of the Appendix.

2.2.7 Constitutive models incorporating machine learning algorithms

There are still some weaknesses for the above traditional nonlinear constitutive models, although intensive studies have been performed and remarkable results have been achieved in the past decades. Most of the physical constitutive models derived from microstructural aspects are expressed in complex equations, but the performances are not always perfect. Moreover, coefficients determined from the tests of both phenomenological models and physical models will affect the performances as well [116]. With the rapid development of machine learning algorithms and computing power, researchers attempt to reflect the behaviors of materials by machine learning algorithms or introducing them into mechanics. In recent years, artificial neural networks (ANN) models are the most popular and widely used machine learning models due to their outstanding performance in terms of “big data” and nonlinear modeling. Therefore, the machine learning-based constitutive models will be reviewed in this section, and the typical model will be recovered and evaluated, so that a clear understanding can be provided of the characteristics of machine learning-based constitutive models. Furthermore, a possible direction for the research of constitutive models can also be provided.

Ghaboussi et al. [118] suggested that the behavior of materials may be reflected by neural networks (NN) as early as 1991. Back-propagation neural network was constructed for the plane stress state under monotonic biaxial loading and compressive uniaxial cycle loading. The NN model is self-operated, and all the parameters responsible for the behavior of materials and experimental data are reflected within a unified environment of a NN. After that, they conducted a series of studies about the NN-based constitutive models, in which the rate-dependent behaviors (as shown in Figure 10) or hysteretic behaviors were considered [117,119]. In their research, the strain, strain rate, and the internal variables were taken as the inputs directly, and stresses were the outputs. The ANN-based model was implemented within the FE method for the results of boundary value problems [120]. The experimental data and simulation data calculated by physical constitutive models were combined, and worked together with the NN model to achieve more accurate results [121]. Their studies not only provided evidence that it is feasible to forecast the behavior of materials by the NN approach, but also established a basic framework for the ANN-based constitutive models.

Figure 10

Rate-dependent NN constitutive model (A _k is equal to e _lm /S _k , where e _lm = ε _lm − δ _lm ε _v /3, S denotes the scale factor, the dots denote the rates, n and n − 1 denote the discrete time steps) [117].

Al-Haik et al. [122] developed an ANN-based constitutive model to predict the stress relaxation of polymer matrix composite. A series of stress relaxation tests were performed at the conditions of constant strain and constant temperature, and 9,000 experimental datasets were used to train the ANN model. The ANN model consisted of three input layers, two hidden layers, and one output layer. To verify the performance of the ANN-based constitutive model, they compared the results of the ANN model and another nonlinear VE constitutive model with the experimental results. The predictions of the ANN model are found to be more accurate than those of the nonlinear VE model over a wider range of stress and temperature conditions, in particular near the T _g. Their study further proved that the ANN model has great potential in constitutive modeling.

Rodriguez et al. [123] predicted the behavior of thermoplastic polymer through a feed-forward artificial NN model. A model with the stress–strain curve as outputs and test datasets as inputs was constructed, and a response surface of stress varying with strain and temperature was modeled as shown in Figure 11. They further compared the performance of ANN-based model with several other constitutive models such as Polynomial model, Neo-Hookean model, Yeoh model, Mooney–Rivlin model, and Ogden model. The proposed ANN-based model achieved the best prediction accuracy among these models, and the predicted maximum percentage error (mean square error) was even less than 1%.

Figure 11

The stress response surface given by the ANN model [123].

By combining the NN model with elasticity theory, researchers also attempted to construct new constitutive models, except for the direct modeling by NN. Jordan et al. [124] decomposed the total logarithmic axial strain into a purely-elastic strain and a VE strain. The elastic deformation was reflected by a temperature-dependent Hooke’s law, while the VE deformation was reflected by the NN model. The relationship between stress with the viscous strain, viscous strain rate, and temperature are identified, without making any prior assumptions on the specific mathematical form. The proposed NN constitutive model and a thermos-elastic-viscoplastic constitutive mode named the Johnsen model [125] were implemented to simulate the deformation of PP at various temperatures and low strain rates, and a comparison was made with experimental results as shown in Figure 12. It was demonstrated that the NN model has the capability to predict the stress–strain response at low strain condition.

Figure 12

Comparison of the NN model and Johnsen model [124].

Li et al. [126] improved the Johnson–Cook model [36] by applying the machine learning algorithm to capture the observed unconventional effect of the strain rate and temperature on the hardening response. A feed-forward NN with three hidden layers was introduced to represent a function that describes the effects of the strain, strain rate, and temperature. The NN was implemented into UMAT, and an accurate NN model was constructed by repeatedly launching full 3D FE simulations during the training. A hardening law of DP800 steel was thus obtained, in which the yield stress is expressed as a function of the equivalent plastic strain, strain rate, and temperature. Although the object of study is metal, this approach is still a typical strategy for combining the NN algorithm with the traditional constitutive model.

Besides, the NN model was not only widely used in the construction of constitutive models, but also in the performance prediction of fiber-reinforced polymeric composite structures [127–129], the damage prediction of composites [130], multiscale modeling of composites [131 or replacing the expensive nonlinear computation of FE model [132].

Tao et al. [133] coupled the commercial FE code Abaqus with the deep neural network model, to acquire the constitutive relationship of the fiber-reinforced composites. The proposed system avoids excessive assumptions in the process of constitutive modeling and satisfies the equilibrium and kinematics equations, so that it has more potential to conform to physics laws. Engineering constants of the fiber‐reinforced composite were acquired, and the progressive damage of structures was also captured. This research shows that the machine learning algorithm can be applied to not only the constitutive modeling of materials, but also the constitutive modeling of structures.

In addition to the above NN algorithm, other machine learning algorithms have also been employed for constructing constitutive models, such as convolutional neural network [134], Gaussian process machine learning [135], and genetic algorithms (GAs) [136]. As far as the authors are concerned, GAs are applied to solve the parameters of constitutive models in most cases. Colak and Cakir [136] proposed a constitutive model to reflect the rate and temperature-dependent stress–strain behavior of polymers. The model was modified from the cooperative-viscoplasticity theory of the VBO model, and GA was adopted to determine the material parameters of the models. The proposed model with parameters obtained from GA was capable of predicting the mechanical behavior and all features below and through T _g of polymers.

The authors built an NN-based constitutive model on the basis of previous research to further evaluate the performance of the machine learning-based constitutive models. A 2D NN model was established, and the model contains one input layer, ten hidden layers, and one output layer. The deformations of PEEK in which strain rate was maintained at 0.001 s⁻¹ were studied, and the experimental data were still taken from the research of Chang et al. [48]. True strain and temperature were the input variables, and true stress was the output variable. The experimental data which obtained from the tests at the temperature of −60℃, −2℃, 23℃, 100℃, and 140℃ were used for training the constitutive model, and the experimental data at the temperature of 60℃ were used for validating the constitutive model. Figure 13 presents the comparison of modeling (surface), training samples (red dots), and validation samples (black dots, the black dots are almost hidden under the surface). It was found that the NN-based constitutive model achieved a good agreement with the experimental results in a wide range.

Figure 13

Modeling of NN-based constitutive model (surface: NN-based model, red dots: samples input, black dots: experimental data to be compared).

In addition, Figure 14 compares the deformation at the temperature of 60℃ predicted with testing separately. Our studies show that the NN-based constitutive model is capable of predicting the deformation of polymers when a suitable model is established with good training accuracy and enough experimental data are provided. in addition, the NN-based constitutive model is suitable for any polymers. Certainly, the performance of a machine learning-based constitutive model greatly depends on the samples, training accuracy, and many other factors, and it is challenging to achieve high accuracy.

Figure 14

Comparison of true stress–strain curves between the prediction of NN-based model with experimental results for PEEK.

Many other researchers were also interested in using the NN algorithms to investigate the nonlinear characteristics of thermoplastic polymers in recent years [137–139]. In addition, the machine learning-based constitutive models in this section are listed in Table A4 of the Appendix. Among the current studies, many benefits are found for such strategies. First, almost no assumptions were made, but the behaviors of materials can be captured by training the NN model by using experimental or computational data. Therefore, such an approach possesses the capability to describe the actual behavior of materials. Second, the NN-based constitutive models may discover unknown deformation mechanisms of polymers by learning the deformation of polymers [116]. Third, the computation may be accelerated by applying the NN algorithm with a reasonable approach.

Up to date, there are still some disadvantages for NN-based constitutive models [116]. In fact, the experiments are usually “expensive,” so it is difficult to obtain a large amount of training data. Due to lack of experimental data, it is not easy to maintain the steady performance of NN-based constitutive models. In addition, the NN-based constitutive models will encounter difficulties to accurately predict the results outside the training sets. For example, when the NN-based model is trained by using experimental data obtained from tensile deformations, such model might encounter difficulties to accurately predict shear deformations. Moreover, the NN-based constitutive model cannot closely relate to the physical mechanism of materials.

Nevertheless, with the development of machine learning technology, these challenges are expected to be overcome in the future. For example, combination of physical constitutive models and machine learning algorithms might be a choice. Therefore, the machine learning algorithms exhibit a great potential in the research of constitutive models for thermoplastic polymers.