*Article* **Thermo-Mechanical Behavior of Poly(ether ether ketone): Experiments and Modeling**

**A. D. Drozdov \* and J. deClaville Christiansen**

Department of Materials and Production, Aalborg University, Fibigerstraede 16, 9220 Aalborg, Denmark; jc@mp.aau.dk

**\*** Correspondence: aleksey@m-tech.aau.dk

**Abstract:** Observations are reported on poly(ether ether ketone) (PEEK) in uniaxial tensile tests, relaxation tests and creep tests with various stresses in a wide interval of temperatures ranging from room temperature to 180 ◦C. Constitutive equations are developed for the thermo–mechanical behavior of PEEK under uniaxial deformation. Adjustable parameters in the governing equations are found by matching the experimental data. Good agreement is demonstrated between the observations and results of numerical simulation. It is shown that the activation energies for the elastoplastic, viscoelastic and viscoelastoplastic responses adopt similar values at temperatures above the glass transition point.

**Keywords:** poly(ether ether ketone); thermo–mechanical response; constitutive modeling

#### **1. Introduction**

Poly(ether ether ketone) (PEEK) is a semicrystalline thermoplastic homopolymer with a linear molecular structure and relatively stiff backbone chains. This polymer belongs to the family of poly(ether ketones) whose ethers functional groups are linked together through aromatic groups. PEEK is a high performance polymer that displays a unique combination of toughness, stiffness, strong abrasion resistance and tribological performance, low moisture absorption, thermo-oxidative stability, chemical and solvent resistance, biocompatibility, flame retardancy, and retention of physical properties at elevated (up to 200 ◦C) temperatures [1]. Due to the excellent balance of mechanical and physical properties, this polymer and composites with PEEK matrices are widely used in aerospace [2,3] and automobile industries [4], energy technologies [5,6] and biomedicine [7,8].

Due to the importance of mechanical properties for application of poly(ether ether ketone) as a load-bearing material, a number of studies have dealt with the experimental investigation of the thermo–mechanical response of PEEK. Results of DMA analysis in the temperature-sweep mode are presented in [9–14]. Observations in tests with various strain rates are reported in [10,15–19] for tensile deformation, in [10,20–22] for compressive deformation, and in [23,24] for biaxial loading. Experimental data in cyclic and fatigue tests are presented in [10,14,25–29]. Observations in the Izod impact tests are given in [13,14,30], whereas those in the Taylor impact tests are provided in [31–33]. Observations in nanoindentation tests are discussed in [34,35]. Experimental data in relaxation tests are reported in [11,36,37], and those in creep tests are presented in [11,38–40].

As PEEK is a semicrystalline polymer, its time- and rate-dependent behavior can be described by conventional models in viscoelasticity and viscoplasticity of semicrystalline polymers, see [41–43], to mention a few. Constitutive equations accounting for the peculiarities in the thermo–mechanical response of PEEK induced by stiffness of its backbone chains were developed in [18,19,22,24,29,44–46].

Previous studies on the thermo–mechanical behavior of PEEK focused on its viscoelastic and viscoplastic responses below the glass transition temperature *T*<sup>c</sup> ≈ 150 ◦C [24,29,44–46]. Above this temperature, only observations in tensile tests with constant strain rates were

**Citation:** Drozdov, A.D.; deClaville Christiansen, J. Thermo-Mechanical Behavior of Poly(ether ether ketone): Experiments and Modeling. *Polymers* **2021**, *13*, 1779. https://doi.org/ 10.3390/polym13111779

Academic Editor: Célio Bruno Pinto Fernandes

Received: 27 April 2021 Accepted: 26 May 2021 Published: 28 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

reported and analyzed in [18,19,22]. The mechanical behavior of PEEK and its micro- and nanocomposites at temperatures exceeding *T*<sup>c</sup> have recently attracted substantial attention due to applications of sulfonated PEEK-based polymers as membrane materials in polymer electrolyte membrane fuel cells and direct methanol fuel cells with the interval of working temperatures up to 180 ◦C [47–49]. The aim of this study is to perform a thorough investigation of the mechanical behavior of PEEK both below and above its glass transition temperature. In particular, we concentrate on the activation energies for the viscoelastic and viscoplastic processes in the sub-*T*<sup>g</sup> and post-*T*<sup>g</sup> intervals. These characteristics allow correlations to be established between mobility of chains and the structure of polymer networks, on the one hand, and physical properties of PEEK membranes in high-temperature electrochemical cells (ionic conductivity, methanol permeability, dielectric permittivity, thermal stability, chemical resistance), on the other [50].

The objective of this paper is three-fold: (i) to analyze experimentally the thermomechanical response of PEEK in uniaxial tensile tests with a constant strain rate, relaxation tests with a constant strain, and creep tests with various stresses in a wide interval of temperatures from room temperature up to 180 ◦C, (ii) to develop constitutive equations for the thermo–elastoplastic, thermo–viscoelastic and thermo–viscoelastoplastic responses of PEEK and to find material constants in these relations by matching the observations, and (iii) to compare activation energies for the elastoplastic behavior (tensile tests with small strains), viscoelastic response (short-term relaxation tests and creep tests in the linear regime of deformation) and viscoelastoplastic behavior (creep tests with relatively large stresses above the glass transition temperature).

Unlike previous studies on modeling the time- and rate-dependent behavior of PEEK subjected to arbitrary 3D deformations with finite strains [19,24,29,44,46], we confine ourselves to the analysis of its thermo-mechanical response under uniaxial deformation with small strains. This allows the number of adjustable coefficients in the constitutive equations to be reduced noticeably (compared with conventional models). As a result, the effect of temperature on the material parameters can be determined with high accuracy. The latter is of primary importance for (i) the design of PEEK implants produced by additive manufacturing, 3D printing and fused deposition modeling (FDM) technology, and (ii) prediction of their microstructure, tribological performance and mechanical properties [51–53].

#### **2. Materials and Methods**

Poly(ether ether ketone) KETRON 1000 PEEK FKM NATUR (density 1.31 g/cm<sup>3</sup> , tensile modulus 4.34 GPa, ultimate tensile strength 110 MPa) was supplied as extruded sheets by Vink Plast ApS (Denmark). Dumbbell specimens for tensile tests (ASTM standard D638) with length in the active zone 50 mm, width 5.1 mm, and thickness 4.5 mm were machined from the sheet. To exclude the effect of stresses developed under preparation, tests were conducted a week after machining samples.

Differential scanning calorimetry (DSC) measurements were carried out by means of STA 449/Netzsch apparatus at the heating and cooling rate of 20 K/min. Specimens with mass of about 10 mg were tested in alumina pans covered by lid under argon atmosphere. The experimental program involves: heating from the initial temperature *T*<sup>i</sup> = 104 up to the final temperature *T*<sup>f</sup> = 400 ◦C, followed by cooling to the initial temperature *T*<sup>i</sup> , and re-heating up to the final temperature *T*<sup>f</sup> . Experimental data are depicted in Figure 1 which shows that the DSC scans for the first and second heating coincide. Figure 1 reveals that the glass transition temperature *T*<sup>g</sup> equals 151 ◦C, the crystallization temperature *T*<sup>c</sup> is 293 ◦C, and the melting temperature *T*<sup>m</sup> equals 339 ◦C. These values are in good agreement with observations in DSC tests on PEEK 90G (*T*<sup>g</sup> = 155, *T*<sup>c</sup> = 317, *T*<sup>m</sup> = 345 ◦C) [14] and PEEK 450G (*T*<sup>g</sup> = 158, *T*<sup>m</sup> = 341 ◦C) [10], as well as with the data in DMA test on PEEK 15G (*T*<sup>g</sup> = 144 ◦C) [11].

**Figure 1.** DSC thermogram of PEEK. Solid line: experimental data in DSC test with the heating and cooling rates of 20 K/min.

Mechanical tests were performed by means of a universal testing machine Instron– 5568 equipped with a thermal chamber and an electro–mechanical sensor (Instron Static 2630113) for control of longitudinal strain in the active zone of samples. Tensile force was measured by a 50 kN load cell. The engineering stress *σ* was calculated as the ratio of axial force to the cross-sectional area of specimens in the undeformed state.

The experimental program included three series of tests at temperatures *T* ranging from room temperature to 180 ◦C. Each test was conducted on a virgin specimen. For each deformation program, tests were repeated three times of different samples to assess repeatability of measurements. The accuracy of measurements is estimated in Supplementary Material (Figures S1–S4) where experimental data are depicted (with their standard deviations) in selected uniaxial tensile tests, relaxation tests, and creep tests together with results of numerical analysis.

The first series involved uniaxial tensile tests with a cross-head speed of 20 mm/min (which corresponded to the strain rate *<sup>e</sup>*˙ <sup>=</sup> 3.1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> s −1 ) up to breakage of specimens. The experimental stress–strain diagrams at temperatures *T* = 20, 80, 120, 130, 140, 150, 160, 170 and 180 ◦C are depicted in Figure 2, where tensile stress *σ* is plotted versus engineering strain *e*. We confine ourselves to the interval 0 ≤ *e* ≤ 0.06 for necking of specimens occur under stretching in the post-yield region at temperatures below the glass transition point *T*g, whereas we focus on the analysis of homogeneous uniaxial deformation with small strains.

**Figure 2.** Stress *σ* versus strain *e*. Symbols: experimental data in tensile tests with strain rate *<sup>e</sup>*˙ <sup>=</sup> 3.1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> s <sup>−</sup><sup>1</sup> at temperatures *T* ◦C. Solid lines: results of simulation.

For each set of data, the maximum stress *σ*max on the stress–strain curve was measured and associated with the yield stress. The effect of temperature on *σ*max is illustrated in Figure 3, where *σ*max is plotted versus *T*. With reference to [54], the data are approximated by the linear equation

$$
\sigma\_{\text{max}} = \sigma\_{\text{max}} \alpha - \sigma\_{\text{max}} \, \_1T\_{\text{\textdegree}} \tag{1}
$$

with the coefficients calculated by the least-squares technique. Figure 3 shows good agreement between the observations and their approximation by Equation (1) with different coefficients below and above the glass transition temperature *T*g.

**Figure 3.** Tensile strength *σ*max versus temperature *T*. Circles: experimental data in tensile tests with strain rate *<sup>e</sup>*˙ <sup>=</sup> 3.1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> s −1 . Solid lines: approximation of the data by Equation (1).

Our findings are in accord with observations on PEEK 450G reported in [10,21], which revealed changes in slope of the dependence *σ*max(*T*) in the interval of temperatures between 135 and 140 ◦C.

The other series of experiments involved tensile relaxation tests with a fixed strain *e*<sup>0</sup> = 0.01. In each test, a specimen was stretched with a cross-head speed of 20 mm/min up to the strain *e*0. Afterwards, the strain was preserved constant, and the tensile stress *σ* was monitored as a function of time *t*. Following the protocol ASTM E-328, a duration of 20 min was chosen for the short-term relaxation tests. Experiments were carried out at temperatures *T* = 20, 80, 120, 130, 140, 150, 160, 170 and 180 ◦C. Selected relaxation curves are reported in Figure 4, where tensile stress *σ* is plotted versus relaxation time *t*rel = *t* − *t*<sup>0</sup> (*t*<sup>0</sup> stands for the time needed to stretch samples up to the strain *e*0). Following common practice, observations are presented by means of the semi-logarithmic plots with log = log10.

**Figure 4.** Stress *σ* versus relaxation time *t*rel. Symbols: experimental data in tensile relaxation tests with strain *e*<sup>0</sup> = 0.01 at temperatures *T* ◦C. Solid lines: results of simulation.

In the last series of experiments, tensile creep tests were performed with various tensile stresses *σ*<sup>0</sup> at various temperatures *T*. In each test, a specimen was stretched with a cross-head speed of 20 mm/min up to the required stress *σ*0. Afterwards, the stress was preserved constant, and an increase in strain *e* was monitored a function of time *t*. Following the protocol ASTM D-2990, a duration of 20 min was chosen for the short-term creep tests.

Two types of creep tests were conducted. Experimental data in these tests are reported in Figures 5 and 6, where tensile strain *e* is plotted versus creep time *t*cr = *t* − *t*<sup>0</sup> (*t*<sup>0</sup> stands for the time needed to reach the stress *σ*<sup>0</sup> under stretching).

In experiments of the first type, tensile stresses *σ*<sup>0</sup> were chosen to be relatively low. Observations in these tests are used to validate our model in linear viscoelasticity. Creep curves in selected tests (with *σ*<sup>0</sup> = 70 MPa at *T* = 20 ◦C, *σ*<sup>0</sup> = 40 MPa at *T* = 120 ◦C and *σ*<sup>0</sup> = 30 MPa at *T* = 150 ◦C) are depicted in Figure 5, and those with *σ*<sup>0</sup> = 10 MPa at *T* = 160, 170 and 180 ◦C are presented in Figure 6.

Tests of the other type were performed with relatively large tensile stresses *σ*<sup>0</sup> to evaluate the viscoplastic flow under creep conditions at temperatures above *T*g. Experimental data in these tests are reported in Figure 6, where creep diagrams are depicted at

temperatures *T* = 160 ◦C (with *σ*<sup>0</sup> = 20, 25, 30 and 35 MPa), *T* = 170 ◦C (with *σ*<sup>0</sup> = 20 MPa) and *T* = 180 ◦C (with *σ*<sup>0</sup> = 30 MPa).

The following conclusions are drawn from Figures 2, 4–6: (i) under tension, stress *σ* decreases monotonically with temperature *T*, (ii) relaxation of stresses is negligible at temperatures below 130 ◦C and becomes noticeable at temperatures above *T*g, (iii) creep flow below the glass transition temperature is weak (an increase in *e* in short-term creep tests does not exceed 0.5%), while this flow becomes pronounced in tests with relatively large stresses above *T*<sup>g</sup> (tensile strain grows by several times).

**Figure 5.** Strain *e* versus creep time *t*cr. Symbols: experimental data in creep tests with various stresses *σ*<sup>0</sup> MPa at temperatures *T* ◦C. Solid lines: results of simulation.

**Figure 6.** Strain *e* versus creep time *t*cr. Circles: experimental data in creep tests with various stresses *σ*<sup>0</sup> MPa at temperatures *T* = 160 ◦C (**A**), *T* = 170 ◦C (**B**) and *T* = 180 ◦C (**C**). Solid lines: results of simulation.

#### **3. Results and Discussion**

We now develop simple constitutive equations for the thermo–mechanical response of PEEK and determine adjustable parameters in these relations by matching the experimental data in Figures 2, 4–6.

#### *3.1. Thermo–Elastoplasticity*

Under tension with a constant strain rate *e*˙, ductile failure (necking) of specimens is observed in the post-yield region. As the failure process is unstable, we study the mechanical behavior of PEEK under tensile deformation up to the points of maximum on the stress–strain diagrams in Figure 2. At all temperatures *T*, tensile stress reaches its ultimate value at strains *e* below 0.05, which corresponds to the duration of stretching of about 16 s. Figure 4 shows that relaxation of stresses during this period does not exceed 7%. This implies that the viscoelastic effects can be disregarded in the analysis of the stress–strain diagrams, and the response of PEEK can be described within the theory of elastoplasticity.

According to this concept, the total strain *e* under uniaxial deformation equals the sum of the elastic, *e*e, and plastic, *e*p, components:

$$
\mathfrak{e} = \mathfrak{e}\_{\mathfrak{e}} + \mathfrak{e}\_{\mathfrak{p}}.\tag{2}
$$

The stress *σ* is connected with the elastic strain *e*<sup>e</sup> by the linear equation

$$
\sigma = \mathbb{E} \mathfrak{e}\_{\mathfrak{e}\prime} \tag{3}
$$

where the Young's modulus *E* is treated as a function of temperature *T*. The plastic strain *e*<sup>p</sup> is connected with the stress *σ* by the flow rule

$$
\epsilon\_\mathbb{P} = A \sinh(B\sigma) \sigma^2,\tag{4}
$$

where *A* is a function of temperature *T*, *B* is a temperature-independent material parameter, and the superscript dot stands for the derivative with respect to time *t*. The multiplier *σ* 2 in Equation (4) provides the simplest version of the Bailey–Norton law (its presence means that the rate of plastic deformation is governed by the stored mechanical energy [41]). The multiplier sinh(*Bσ*) is introduced to avoid the use of a yield surface (plastic deformation is presumed to occur at any stress, but its rate is negligible at small stresses due to the properties of the hyperbolic sinus). Another explanation for this term is based on the Eyring theory of thermally activated processes, according to which *B* is proportional to the activation volume for cooperative motion of polymer chains [55].

Equations (2)–(4) provide constitutive relations in thermo–elastoplasticity of PEEK under uniaxial deformation. These equations involve two functions of temperature, *E*(*T*) and *A*(*T*), and a constant *B* to be found by approximation of experimental data in Figure 2.

We begin with matching observations in tensile tests at room temperature, calculate *E* by fitting the data at 0 ≤ *e* ≤ 0.01, and determine *A* and *B* from the best-fit condition for the entire stress–strain diagram. Then, the coefficient *B* = 0.03 MPa−<sup>1</sup> is fixed, and the stress–strain curves at elevated temperatures are matched by means of two parameters, *E* and *A*, only. Each set of observations is approximated separately. Figure 2 shows an acceptable agreement between the data and results of simulation.

To evaluate the activation volume *V*<sup>a</sup> associated with the coefficient *B* and to compare its values with results of other studies, the theory of thermally activated processes is applied. According to this concept,

$$B = \frac{V\_\mathrm{a}\sqrt{3}}{k\_\mathrm{B}T\_0},\tag{5}$$

where *k*<sup>B</sup> is the Boltzmann constant, *T*<sup>0</sup> = 293 K stands for room temperature, and the coefficient <sup>√</sup> 3 appears due to transformation of tensile stress into the equivalent shear stress. For *B* found by fitting observations in Figure 2, it follows from Equation (5) that *<sup>V</sup>*<sup>a</sup> <sup>=</sup> 7.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> nm<sup>3</sup> . This value is substantially (by two to three orders of magnitude) lower than *V*<sup>a</sup> = 1 nm<sup>3</sup> [16], *V*<sup>a</sup> = 1 to 7 nm<sup>3</sup> [21], *V*<sup>a</sup> = 3.4 nm<sup>3</sup> [35], and *V*<sup>a</sup> = 12.6 nm<sup>3</sup> [30]. This difference can be attributed to the fact that the activation volumes were calculated in the above works as measures of sensitivity of the yield stress to changes in the strain rate [55]. This explains also large (by an order of magnitude) deviations between the reported values of *V*a.

The situation changes drastically if we associate *V*<sup>a</sup> with volumes of holes measured by means of the positron annihilation lifetime spectroscopy (PALS) and diffusivity of gases (these two methods lead to similar results [56]). PALS measurements of free volume imply that *<sup>V</sup>*<sup>a</sup> <sup>=</sup> 7.15 <sup>×</sup> <sup>10</sup>−<sup>2</sup> nm<sup>3</sup> for PEEK specimens [57], and this parameter varies between 7.2 <sup>×</sup> <sup>10</sup>−<sup>2</sup> to 8.4 <sup>×</sup> <sup>10</sup>−<sup>2</sup> nm<sup>3</sup> depending on degree of crystallinity [58]. Both estimates are in good accord with the value of *V*<sup>a</sup> obtained in our analysis of observations.

The effects of temperature *T* on the elastic modulus *E* and the rate of plastic flow *A* are illustrated in Figures 7 and 8. By analogy with Equation (1), the data in Figure 7 are approximated by the linear function

$$E = E\_0 - E\_1 T \tag{6}$$

with the coefficients accepting different values below and above the glass transition temperature *T*g. Figure 7A shows that the modulus *E* remains practically constant below the glass transition temperature and decreases strongly above *T*g. This behavior resembles that observed for the shear storage modulus [10] and tensile storage modulus [9,11] in DMA tests on PEEK.

**Figure 7.** The Young's modulus *E* versus temperature *T*. Symbols: (**A**)—treatment of observations in tensile tests (◦) and creep tests (•). (**B**)—treatment of observations in relaxation tests. Solid lines: approximation of the data by Equation (6).

The coefficient *A* is plotted versus reciprocal absolute temperature in Figure 8. The data are approximated by the Arrhenius dependence

$$A = A\_0 \exp\left(-\frac{E\_\mathrm{a}}{RT}\right),\tag{7}$$

where *A*<sup>0</sup> is a pre-factor, *R* is the universal gas constant, and *E*<sup>a</sup> stands for the activation energy. Figure 8 shows that the observations are correctly described by Equation (7) when different activation energies are used below and above the glass transition temperature *T*g.

**Figure 8.** Coefficient *A* versus temperature *T*. Circles: treatment of observations in tensile tests. Solid lines: approximation of the data by Equation (7) with *E*<sup>a</sup> = 16.7 kJ/mol (low temperatures) and *E*a = 133.3 kJ/mol (high temperatures).

#### *3.2. Thermo–Viscoelasticity*

The experimental data in tensile relaxation tests with a small strain *e*<sup>0</sup> = 0.01 at various temperatures *T* are described by the constitutive equations in linear viscoelasticity of semicrystalline polymers [59]. A polymer is thought of as a network with two types of chains: permanent (whose ends are bridged by covalent cross-links) and temporary (able to separate from their junctions and to merge with the network at random times being driven by thermal fluctuations). The heterogeneous network is composed of meso-domains with various activation energies for rearrangement. The rate Γ for separation of temporary chains from their junctions in meso-domains with a dimensionless activation energy *u* (normalized by *k*B*T*0) is governed by the Eyring equation

$$
\Gamma = \gamma \exp(-\mu),
\tag{8}
$$

where *γ* is a pre-factor. A quasi-Gaussian expression is adopted for the distribution function *f*(*u*) of meso-domains with various activation energies *u*,

$$f(u) = f\_0 \exp\left(-\frac{u^2}{2\Sigma}\right) \qquad (u \ge 0). \tag{9}$$

The dimensionless parameter Σ characterizes inhomogeneity of an ensemble of mesodomains. The coefficient *f*<sup>0</sup> is determined by the normalization condition

$$\int\_0^\infty f(u)\mathbf{d}u = 1.\tag{10}$$

Under uniaxial tension with an arbitrary deformation program *e*(*t*), tensile stress *σ*(*t*) obeys the constitutive equation

$$\sigma(t) = E\left[\varepsilon(t) - \kappa \int\_0^\infty \Gamma(u) f(u) \mathrm{d}u \int\_0^t \exp\left(-\Gamma(u)(t-\tau)\right) \varepsilon(\tau) \mathrm{d}\tau\right],\tag{11}$$

where *E* stands for the Young's modulus, and *κ* denotes the ratio of the number of temporary chains to the total number of chains per unit volume in the initial state. Unlike conventional (the Maxwell–Wiechert type) models for the linear viscoelastic response of polymers (involving a large number of material constants), Equation (11) is entirely determined by four parameters, *E*, *κ*, *γ* and Σ. With reference to [60], we suppose that Σ is independent of temperature, *κ* increases linearly with temperature and reaches its ultimate value *κ* = 1 at relatively high temperatures *T* > *T*g, and *γ* obeys the Arrhenius law. Equation (11) implies that, in tensile relaxation tests with a fixed strain *e*0, the stress *σ* decreases with relaxation time *t*rel following the pattern

$$\sigma(t\_{\rm rel}) = \sigma\_0 \left\{ 1 - \kappa \int\_0^\infty f(u) \left[ 1 - \exp\left( -\gamma \exp(-u) t\_{\rm rel} \right) \right] du \right\},\tag{12}$$

where

$$
\sigma\_0 = E \epsilon\_0 \tag{13}
$$

stands for the stress at the beginning of the relaxation process.

Each set of data in Figure 4 is matched separately by means with Equation (12) with three parameters *σ*0, *κ* and *γ* (since Σ is independent of temperature, we find its value Σ = 7.0 by fitting observations at *T* = 160 ◦C and use it at all temperatures *T* under consideration). The best-fit coefficient *γ* is determined by means of the nonlinear regression algorithm, while *σ*<sup>0</sup> and *κ* are determined by the least-squares technique. Given *σ*0, the modulus *E* is calculated from Equation (13) and plotted versus temperature *T* in Figure 7B. This figure demonstrates that the data are described adequately by Equation (6) with different coefficients *E*<sup>0</sup> and *E*<sup>1</sup> below and above the glass transition temperature *T*g.

The coefficients *E*<sup>0</sup> and *E*<sup>1</sup> found by approximation of experimental data in Figure 7A,B above *T*<sup>g</sup> coincide practically: the difference is less than 7% for *E*<sup>0</sup> and 5% for *E*1. The discrepancies between the coefficients *E*<sup>0</sup> and *E*<sup>1</sup> calculated by matching observations in tensile and relaxation tests below *T*<sup>g</sup> are higher (but do not exceed 14%). They may be explained by local variations in thicknesses of specimens machined from an extruded sheet.

The parameter *κ* is plotted versus temperature in Figure 9A. This figure demonstrates that *κ* increases monotonically with *T* (the growth of intensity of thermal fluctuations results in transformation of some permanent chains into transient chains). The influence of temperature *T* on *κ* is described by the linear equation

$$
\kappa = \kappa\_0 + \kappa\_1 T\_\prime \tag{14}
$$

where *κ*<sup>0</sup> and *κ*<sup>1</sup> are calculated by the least-squares method. Figure 9A shows that Equation (14) describes correctly the function *κ*(*T*) when different coefficients *κ*0, *κ*<sup>1</sup> are used below and above the glass transition temperature *T*g.

The effect of temperature *T* on the rate of relaxation *γ* is illustrated in Figure 9B. The data (in the region of temperatures *T* ≥ 120 ◦C) are approximated by the Arrhenius equation

$$\gamma = \gamma\_0 \exp\left(-\frac{E\_\mathrm{a}}{RT}\right),\tag{15}$$

where *γ*<sup>0</sup> stands for a pre-factor, and *E*<sup>a</sup> denotes the activation energy. Comparison of Figures 8 and 9B shows that the activation energies found by matching observations in tensile tests and relaxation tests in the high temperature region adopt similar values (deviations between *E*<sup>a</sup> = 133.3 kJ/mol in Figure 8 and *E*<sup>a</sup> = 113.4 kJ/mol in Figure 9B do not exceed 15 %).

**Figure 9.** Parameters *κ* and *γ* versus temperature *T*. Circles: treatment of observations in relaxation tests. Solid lines: (**A**)—approximation of the data by Equation (14). (**B**)—approximation of the data by Equation (15) with *E*a = 113.4 kJ/mol.

It is worth noting that these activation energies differ pronouncedly from the activation energies for *α*-relaxation reported in previous studies on the time-dependent response of PEEK (*E*<sup>a</sup> = 377 [38], *E*<sup>a</sup> = 494 [37] *E*<sup>a</sup> = 582 [61], *E*<sup>a</sup> = 810 to 1000 [9], *E*<sup>a</sup> = 1094 [39], and *E*<sup>a</sup> = 1100 kJ/mol [62]). To explain this difference, experimental data in relaxation tests are treated by the conventional method: all relaxation curves at elevated temperatures *T* are shifted horizontally (along the time-axis) to construct a master-curve at room temperature *T*0. Figure 10A confirms that a smooth master-curve is formed by means of this technique. The shift factor *a* is plotted versus temperature *T* in Figure 10B. The data are approximated by the Arrhenius equation

$$
\log a = \log a\_0 + \frac{E\_\mathrm{a}}{RT} \,\mathrm{}\,\tag{16}
$$

where *a*<sup>0</sup> and *E*<sup>a</sup> are calculated by the least-squares technique.

**Figure 10.** (**A**)—Stress *σ* versus relaxation time *t*rel. Symbols: master-curve at room temperature constructed by shift of experimental data in relaxation tests. (**B**)—Shift factor *a* versus temperature *T*. Circles: treatment of observations in relaxation tests. Solid lines: approximation of the data by Equation (16) with *E*<sup>a</sup> = 159.4 kJ/mol (low temperatures) and *E*a = 531.3 kJ/mol (high temperatures).

Figure 10B shows good agreement between the data and predictions of Equation (16) with different coefficients below and above the glass transition temperature *T*g. The

activation energy in the low-temperature region *E*<sup>a</sup> = 159.4 kJ/mol is in accord with the activation energy *E*<sup>a</sup> = 177 kJ/mol for the *β*-relaxation in PEEK [9], whereas the activation energy in the high-temperature region *E*<sup>a</sup> = 531.3 kJ/mol is close to the activation energies *E*<sup>a</sup> = 494 kJ/mol [37] and *E*<sup>a</sup> = 582 kJ/mol [61] for the *α*-relaxation. These results demonstrate that the method based on the time-temperature superposition principle (shifts of observations in relaxation and dynamic mechanical tests) overly estimates the activation energies since this approach disregards evolution of *E* and *κ* with temperature *T*.

To confirm validity of our model for the linear viscoelastic behavior of PEEK, we apply Equation (11) to describe the time-dependent response of specimens in tensile creep tests and compare results of simulation with experimental data depicted in Figures 5 and 6. Resolving Equation (11) with respect to *e*, we find that the growth of tensile strain with time in a creep test with a constant stress *σ*<sup>0</sup> is governed by the equation

$$
\epsilon(t\_{\rm cr}) = \epsilon\_0 + \kappa \int\_0^\infty f(u) s(t\_{\rm cr}, u) \mathrm{d}u. \tag{17}
$$

Here *t*cr = *t* − *t*0, where *t*<sup>0</sup> is the moment when tensile stress reaches the required value *σ*0,

$$
\epsilon\_0 = \frac{\sigma\_0}{E} \tag{18}
$$

is the strain at the beginning of the creep process, and the function *s*(*t*cr, *u*) obeys the differential equation

$$\frac{\partial s}{\partial t\_{\rm cr}}(t\_{\rm cr}, u) = \gamma \exp(-u) \left[ \varepsilon(t\_{\rm cr}) - s(t\_{\rm cr}, u) \right], \qquad s(0, u) = 0. \tag{19}$$

For each temperature *T* and stress *σ*0, Equations (17)–(19) are integrated numerically by the Runge–Kutta method with the material constants *γ*, *κ* and Σ determined by matching observations in relaxation tests (Figure 9). To ensure good agreement between the observations in Figures 5 and 6 and results of simulation, the coefficient *E* is treated as an adjustable parameter. Its best-fit values at various temperatures *T* are reported in Figure 7A. This figure shows that the Young's moduli found by matching experimental data in tensile tests and creep tests coincide practically. Comparison of experimental data with results of simulation (Figures 5 and 6) confirms the ability of the model (17)–(19) to predict the response of PEEK in short-term creep tests with small strains when material parameters are determined in relaxation tests.

#### *3.3. Thermo–Viscoelastoplasticity*

The time-dependent response of PEEK in creep tests with relatively high stresses *σ*<sup>0</sup> (beyond the interval of linear viscoelasticity) is described within the concept of viscoelastoplasticity [63]. In accord with Equation (2), the total strain *e* is split into the sum of the viscoelastic strain, *e*ve, and plastic strain, *e*p,

$$
\mathfrak{e} = \mathfrak{e}\_{\text{ve}} + \mathfrak{e}\_{\text{p}}.\tag{20}
$$

The viscoelastic strain *e*ve is connected with tensile stress *σ* by Equation (11),

$$\sigma(t) = E\left[\mathfrak{e}\_{\rm ve}(t) - \mathfrak{x} \int\_0^\infty \Gamma(u) f(u) \mathrm{d}u \int\_0^t \exp\left(-\Gamma(u)(t-\tau)\right) \mathfrak{e}\_{\rm ve}(\tau) \mathrm{d}\tau\right]. \tag{21}$$

The evolution of the plastic strain *e*<sup>p</sup> with time *t* is governed by an analog of Equation (4),

$$
\dot{\epsilon}\_{\rm p} = \bar{A}\sinh(\epsilon)(\sigma - \sigma\_{\rm b})^{m}, \qquad \epsilon\_{\rm p}(0) = 0,\tag{22}
$$

where *A*¯ stands for a pre-factor, the term sinh(*e*) is introduced to avoid the growth of plastic strain at small strains (far below the yield point), and *m* is an exponent in the Bailey–Norton law. The back-stress *σ*<sup>b</sup> accounting for strain hardening under plastic flow [64] reads

$$
\sigma\_{\mathsf{b}} = E\_{\mathsf{b}} \epsilon\_{\mathsf{p}\prime} \tag{23}
$$

where *E*<sup>b</sup> is an analog of the elastic modulus.

Equations (20)–(23) together with Equations (8) and (9) provide a constitutive model in viscoelastoplasticity of semicrystalline polymers under uniaxial deformation. These relations involve four parameters, *E*, *κ*, *γ* and Σ, that characterize the linear viscoelastic response, and three extra coefficients, *A*¯, *E*<sup>b</sup> and *m*. To reduce the number of parameters, we set *m* = 7 (a typical value of the Bailey–Norton exponent) in matching experimental creep curves in Figure 6. The other two quantities, *A*¯ and *E*b, are considered as functions of temperature *T* only.

Under tensile creep with a fixed stress *σ*0, Equation (21) takes a form similar to Equation (17),

$$
\epsilon\_{\rm ve}(t\_{\rm cr}) = \epsilon\_0 + \kappa \int\_0^\infty f(u) s(t\_{\rm cr}, u) \mathrm{d}u,\tag{24}
$$

while Equation (22) yields

$$\frac{d\epsilon\_{\rm p}}{dt\_{\rm cr}} = \bar{A}\sinh(\varepsilon)(\sigma\_0 - \sigma\_{\rm b})^m, \qquad \epsilon\_{\rm p}(0) = 0. \tag{25}$$

The following procedure is applied to fit observations in Figure 6 at each temperature under consideration. First, the creep diagram with the highest stress *σ*<sup>0</sup> is approximated with the help of two parameters, *A*¯ and *E*<sup>b</sup> (the quantities *E*, *κ*, *γ* and Σ are found by matching experimental data in the corresponding relaxation test). Afterwards, the quantities *A*¯, *E*b, *m*, *κ*, *γ*, Σ are fixed, and each remaining creep curve is fitted with the only coefficient *E*. We treat *E* as an adjustable parameter to account for deviations in thicknesses of specimens machined from an extruded sheet. Discrepancies between the best-fit values of *E* determined by matching creep curves with various *σ*<sup>0</sup> have the same order of magnitude as those between the data in tensile and creep tests in Figure 7 (they do not exceed 12 % at all temperatures).

The effects of temperature on coefficients *A*¯ and *E*<sup>b</sup> are illustrated in Figure 11.

**Figure 11.** Parameters *A*¯ and *E*<sup>b</sup> versus temperature *T*. Circles: treatment of observations in creep tests. Solid lines: (**A**)—approximation of the data by Equation (26) with *E*a = 147.6 kJ/mol. (**B**)—approximation of the data by Equation (27).

The growth of the kinetic parameter *A*¯ with *T* is described by the Arrhenius dependence,

$$
\bar{A} = \bar{A}\_0 \exp\left(-\frac{E\_\mathrm{a}}{RT}\right). \tag{26}
$$

Figure 11A shows an acceptable agreement between the data and their predictions by Equation (26) with the activation energy *E*<sup>a</sup> = 147.6 kJ/mol. The latter value is in agreement with the activation energy *E*<sup>a</sup> = 133.3 kJ/mol determined by matching observations in Figure 8 (the difference is less than 10 %).

The decay in *E*<sup>b</sup> with temperature is described by the equation analogous to Equation (6),

$$E\_{\rm b} = E\_{\rm b0} - E\_{\rm b1}T\_{\rm \prime} \tag{27}$$

where *E*b0 and *E*b1 are found by the least-squares method. Figure 11B demonstrates good agreement between the data and results of simulation. A similarity should be stressed between the effect of temperature on the Young's modulus *E* and the modulus *E*<sup>b</sup> that characterize the back-stress. Above the glass transition temperature *T*g, the dimensionless ratios *<sup>E</sup>*1/*E*<sup>0</sup> adopt the values 5.09 <sup>×</sup> <sup>10</sup>−<sup>3</sup> (Figure 7A) and 5.01 <sup>×</sup> <sup>10</sup>−<sup>3</sup> (Figure 7B), whereas the ratio *<sup>E</sup>*b1/*E*b0 equals 4.62 <sup>×</sup> <sup>10</sup>−<sup>3</sup> (Figure 11B), which confirms that changes in *E* and *E*<sup>b</sup> with temperature are governed by the same physical mechanism.

#### **4. Conclusions**

Observations are reported on poly(ether ether ketone) KETRON 1000 PEEK in DSC tests (with a constant rate of 20 K/min under heating and cooling), as well as in uniaxial tensile tests with constant strain rate *<sup>e</sup>*˙ <sup>=</sup> 3.1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> s −1 , relaxation tests with a constant strain *e*<sup>0</sup> = 0.01, and creep tests with various stresses *σ*<sup>0</sup> (ranging from 10 to 70 MPa) at temperatures ranging from room temperature to 180 ◦C.

Constitutive equations are developed for the thermo–elastoplastic, thermo–viscoelastic and thermo–viscoelastoplastic responses of PEEK under uniaxial deformation with small strains. An advantage of these relations is that they involve relatively small numbers of adjustable parameters (three for the elastoplastic behavior, four for the linear viscoelastic response, and an extra three for the viscoelastoplastic flow). Good agreement is demonstrated between the experimental data and results of simulation.

Analysis of the effect of temperature on the rate of elastoplastic deformation *A*, the relaxation rate *γ*, and the rate of viscoplastic flow *A*¯ reveals that the growth of these quantities with *T* above the glass transition temperature *T*<sup>g</sup> obeys the Arrhenius law with similar activation energies *E*<sup>a</sup> (in the range between 113 and 148 kJ/mol).

The study of the elastoplastic and viscoelastoplastic responses of PEEK shows that (i) the activation volume for plastic deformation coincides with the free volume found in PALS and diffusion tests and (ii) the elastic moduli *E* (stress) and *E*<sup>b</sup> (back-stress) decrease similarly with temperature.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/polym13111779/s1, Figure S1: Stress *σ* versus strain *e*. Symbols: experimental data in tensile tests at temperatures *T* = 20 (◦) and 170 (•) ◦C. Bars stand for the standard deviations. Solid lines: results of simulation, Figure S2: Stress *σ*<sup>0</sup> versus relaxation time *t*rel. Symbols: experimental data in tensile relaxation tests with strain *e*<sup>0</sup> = 0:01 at temperatures *T* = 20 (◦) and 170 (•) ◦C. Bars stand for the standard deviations. Solid lines: results of simulation, Figure S3: Strain *σ* versus creep time *t*cr. Circles: experimental data in creep test with stress *σ*<sup>0</sup> = 70:0 MPa at temperature *T* = 20 ◦C. Bars stand for the standard deviations. Solid line: results of simulation, Figure S4: Strain *e* versus creep time tcr. Symbols: experimental data in creep tests with stresses *σ*<sup>0</sup> = 10 (◦) and 20 (•) MPa at temperature *T* = 170 ◦C. Bars stand for the standard deviations. Solid lines: results of simulation.

**Author Contributions:** A.D.D.: Conceptualization, methodology, software, formal analysis, data curation, writing—original draft preparation, writing—review and editing. J.d.C.: Conceptualization, methodology, investigation, writing—original draft preparation, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported Innovationsfonden (Innovation Fund Denmark).

**Institutional Review Board Statement:** not applicable

**Informed Consent Statement:** not applicable

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** Financial support by Innovationsfonden (Innovation Fund Denmark, project 9091-00010B) is gratefully acknowledged.

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

#### **References**

