Study of the Thermomechanics of the Additive Manufacturing Process of Biocompatible Products Subject to the Viscoelastic Behavior of the Functional Material Polyetheretherketone

Abstract

This study considers the problem of numerical modeling of the PEEK product’s 3D printing using the FDM technology. The aim of the study is to verify the adequacy of the use of a thermoviscoelastic model for numerical computations of the PEEK deposition process and to develop an algorithm for calculating this process. The Prony model is used to describe the thermoviscoelastic behavior of the material under study; the temperature-time shift is described by the Williams–Landel–Ferry function (WLF). To obtain the values of the material constants of the relaxation function, first, we used data from other authors; however, after their substitution into the numerical simulation, it was not possible to obtain results close to the full-scale experiment. Therefore, realized our own DMA experiment. The algorithm was developed and implemented in the ANSYS package to calculate non-stationary temperature fields and the stress–strain state of the structure during its layer-by-layer deposition. To solve these problems, the technology of “killing” and subsequent “aliving” of the PEEK material, implemented in the ANSYS package, is used. The numerical algorithm is verified with the results of an experiment on printing samples from PEEK. A good consistency of the calculated data with the experiment is shown.

1. Introduction

Polyetheretherketone (PEEK) is a thermoplastic semi-crystalline polymer that is actively used in modern industries due to a unique combination of a number of operational properties [1], thermo, heat, and fire resistance, low hygroscopicity, radiation and chemical resistance, good mechanical and dielectric properties, and biocompatibility [2,3,4].

It is known from [5,6] that the crystallization temperature of PEEK, depending on the manufacturer, varies from 287 to 311 °C, the melting temperature varies from 339 to 347 °C, the decomposition temperature varies from 567 to 589 °C, and the glass transition temperature varies from 140 to 148 °C, the upper limits of the operating temperatures being from 230 to 290 °C, respectively.

It is known from [7] that the PEEK density is a function of the degree of crystallinity and can be assumed constant for a material in the amorphous ρ_a = 1262.6 Kg/m³ and crystallized state ρ_k = 1400.6 Kg/m³.

Similarly to [7], the coefficient of thermal expansion can be considered a constant value at the temperature intervals of the glassy α_c = 100 × 10⁻⁶ 1/°C and highly elastic states α_b = 670 × 10⁻⁶ 1/°C.

The dependences of thermophysical characteristics of PEEK on temperature are known from [8,9] and represent non-monotonic curves over the studied temperature range from 25 to 330 °C. In the range from 25 to 100 °C, a decrease in the specific heat capacity from C_p = 1242 J/(kg K) to 1.181 J/(kg K) and a subsequent return to the previous level is observed for the curves. Similarly, for thermal conductivity: from λ = 0.268 W/(m K) to 0.157 W/(m K). In the temperature range from 100 to 300 °C, an almost linear increase in the described characteristics is observed from C_p = 1282 J/(kg K) to 2412 J/(kg K) and from λ = 0.230 W/(m K) to 0.264 W/(m K), respectively. When the crystallization temperature is reached, a sharp increase in both characteristics is observed up to the limit values of C_p = 3248 J/(kg K) and λ = 0.357 W/(m K) with their further decrease as the melting temperature is reached. At the same time, when solving a number of problems (for example, calculating the stress–strain state of the structure [7]), the mentioned thermophysical properties of the material can be assumed constant and equal to C_p = 1300 J/ (kg K) and λ = 0.25 W/(m K).

There are a number of comparative calculations on the mechanical behavior of PEEK in extreme operating conditions (with the transition to plastic) [10]. To describe the thermal deformation of the PEEK material and calculate the stress–strain state of the structure, the Johnson–Cook model is mainly used. This model allows taking into account plastic deformation, strain rate, viscous effects, and the thermal softening of the material [11,12,13,14,15]. However, this model has some drawbacks. It is convenient to use in software packages with an explicit scheme for calculating non-stationary problems. We had an extremely poor convergence of the solution when trying to implement it in ANSYS Mechanical APDL. In addition, the validity of using this model to describe polymers remains an open question. Moreover, it should be noted that in the process of layer-by-layer growth, the material works under normal conditions, not associated with large active deformations, but in a wide temperature range. In this case, the main hypothesis of the formation of residual deflections [7] is the assumption of the prevailing contribution of the incompatibility of deformations in individual layers due to the non-synchronous solidification of the thermoviscoelastic material. A similar mechanism is used in the production of prestressed glass. For this reason, a viscoelastic model was chosen.

The Prony model assumes a preliminary determination of the values of the material constants of the relaxation function. As traditionally, the data from the tests for uniaxial tension–compression, not for shear, are used to obtain experimental parameters, the shear constants are calculated through relations linking them to constants from the tension–compression experiment, which are found from the experiment to determine the complex modulus (Section 2.2 and Section 2.2.1).

Numerical values of the temperature dependence function of the complex relaxation modulus E can be calculated on the basis of experimental data from [7,16,17,18] (Section 2.2.3). Experimental graphs of the amplitude of the complex relaxation E modulus at sample cooling rates of 1, 10, and 35 °C/s known from [7] demonstrate a sharp increase in the modulus from 5 × 10⁷ Pa to 3.2 × 10⁹ Pa at a temperature of 173 °C and from 5 × 10⁷ Pa to 2.7 × 10⁹ Pa at a temperature of 143 °C, respectively. The model for predicting the reaction of a material during technological cooling, described in these sources, includes temperature evolution, crystallization kinetics, and a viscoelastic model for predicting thermomechanical properties. The kinetics of crystallization in semi-crystalline thermoplastic composites is considered in connection with the influence of the degree of crystallinity on the mechanical properties and possible contribution of volumetric shrinkage deformation, and the model of the nucleation and growth of crystals of Velisaris and Seferis is used to describe it [16,17]. To predict dynamic modules in [7], a modified form of the standard linear viscoelastic solid model was used. This model was expanded to account for the effect of crystallinity on the behavior of semi-crystalline thermoplastic matrices. Experimental data describing the process of crystallization of PEEK are known from [7,16,17]. Study [7] presents the results of modeling the dependence of the material crystallinity degree on temperature for cooling rates of 1, 10, and 35 °C/s. A sharp decrease in the degree of crystallinity of the material is observed from 0.27 to 0 at 280 °C, from 0.24 to 0 in the temperature range from 200 to 260 °C, and from 0.025 to 0 in the temperature range from 200 to 250 °C for the corresponding cooling rates. The experimental curves of the degree of material crystallinity for cooling rates of 9.4 °C/min, 19.2 °C/min, 37.1 °C/min, and 55.8 °C/min are given in [16,17]. A decrease in the material crystallinity degree is observed from 0.32 to 0 in 250 s, from 0.3 to 0 in 100 s, from 0.28 to 0 in 70 s, and from 0.26 to 0 in 25 s for the corresponding cooling rates. Experimental curves describing the dependence of the material crystallinity degree on time at set temperatures of 307, 310, 312, and 315 °C are also given in [16,17]. The material crystallinity degree reaches values of 0.3 per 1000 s, 0.34 per 1600 s, and 0.33 per 2100 and 3200 s, respectively, at the corresponding temperatures. The dependence of the elastic modulus of the material on the crystallinity degree of the material of the sample is also known from [19]. There is an increase in the elastic modulus of the material from 2.8 to 5.5 GPa with an increase in the material crystallinity degree from 0 to 35%, respectively.

The numerical modeling of the process of layer-by-layer deposition of polymer materials is considered in [20,21,22,23]. The existing solutions cover all stages of the extrusion deposition process of the material: the feeding, melting, extrusion, and deposition of the substance. Models of the material feeding process allow us to calculate the optimal feed rate of the solid filament and the pushing force, as well as to take into account and minimize the effect of slipping of the filament. Models of the material melting process describe the viscous behavior of the melt, the dependence of viscosity on temperature, internal heat exchange, and heat exchange with the surrounding medium, as well as the effect of the nozzle angle and the material feed rate on pressure changes. Models of the material extrusion process allow us to calculate the expansion and convective cooling of the material when exiting the nozzle, and models of the material deposition process—the spreading, cooling, and adhesion of the material. The results of these studies indicate that numerical modeling is in good agreement with experimental values when predicting the properties of materials.

Furthermore, there are a number of studies in which, in a similar way, using the technology of the birth and death of elements, the processes of melting of wire material are simulated. The model takes into account the thermal and deformation processes observed in the sample, and the data of the computational experiment are in good agreement with the full-scale experiment [24,25].

In [26,27,28], the influence of printing parameters on the change of the finished product shape is analyzed, and their optimal values for PEEK are given.

2. Materials and Methods

2.1. Mathematical Formulation

Taking into account small deformations and negligibly small dissipative heat release, it is possible to separate the boundary value problem of unsteady thermal conductivity and the boundary value problem of thermomechanics on the stress–strain state (SSS), which are unrelated in this formulation. The technology of “killing” and subsequent “aliving” (Elements Birth and Death technology in ANSYS) of a part of the material that was initially absent in the model and appearing during the deposition process is used to solve them. At the same time, the domain occupied by the finished product is considered as the computational one. The continuous deposition of the material is carried out discretely; at each calculation sub-stage corresponding to the “aliving” of the next subdomain of the “killed” elements, the boundary value problem of thermal conductivity and thermomechanics is solved, and the result of solving the previous sub-stage serves as the initial conditions for the next one (see Section 2.3.1 and Section 2.3.2 for a detailed illustration of the simulation of the material deposition process). Modeled geometry is presented in Figure 1.

At the $k$-th sub-stage of the solution, the formulation of the boundary value problem of unsteady thermal conductivity for determining temperature fields $T\left(x,t\right)$ in a domain $V_k$ with a boundary $S_{1,k}{{\displaystyle \cup }}^{ }{S}_2$, taking into account the accepted hypotheses, includes [29]:

heat conductivity equation:

(1)$\rho \left(x\right)c\left(x,T\right)\frac{\partial T}{\partial t}=\mathrm{div}\left(\mathsf{\lambda }\left(x,T\right)\mathrm{grad}\left(T\right)\right)+\rho \left(x\right)\dot{q}\left(x,t\right), x\in V_k,$

where $c\left(x,T\right)$, $\mathsf{\lambda }\left(x,T\right)$, and $\rho \left(x\right)$ are the heat capacity, heat conductivity, and density of the material, respectively. $\dot{q}(x,t)$ is the specific power of an external heat source, which means the energy that is taken by the extruder heater for heating PEEK up to 430 °C;

boundary conditions:

(2)$-\mathsf{\lambda }\left(x,T\right)\mathrm{grad}\left(T\right)\cdot n=h\left(T\right)\cdot \left(T-T_c\left(t\right)\right),x\in S_{1,k},$

(3)$T\left(x,t\right)=T_o\left(x\right), x\in S_2,$

where the first equation describes convective heat transfer, and the second one describes platform heating; $h\left(T\right)$ is the heat transfer coefficient, $T_c\left(t\right)$ is the ambient temperature, $n$ is the normal vector to the boundary $S_{1,k}$ of objects, and $T_o\left(x\right)$ is the platform heating temperature;

initial conditions:

(4)$T\left(x,t_{0,k}\right)=T_{k-1}\left(x\right), x\in V_k,$

where $T\left(x,t_{0,k}\right)$ is the initial temperature distribution for the $k$-th sub-stage and $T_{k-1}\left(x\right)$ is the temperature determined at the end of the previous one.

These relations take into account the condition that the research domain $V_k=V_k^{liv}{{\displaystyle \cup }}^{ }{V}_k^{kil}$ remains unchanged throughout the sub-stage. Here, $V_k^{liv}$ and $V_k^{kil}$ indicate zones occupied by “alive” and “killed” elements, respectively. At the same time, the thermophysical properties of the material in the zone of “killed” elements are subject to degradation:

$$\begin{gathered} c\left(x\right), x\in V_k^{kil}<<c\left(x,T\right), x\in V_k^{liv}, \rho \left(x\right), x\in V_k^{kil}<<\rho \left(x,T\right), x\in V_k^{liv}, \\ \mathsf{\lambda }\left(x\right), x\in V_k^{kil}>>\mathsf{\lambda }\left(x,T\right), x\in V_k^{liv}. \end{gathered}$$

The unrelated quasi-static boundary value problem of solid mechanics taking into account the insignificance of the contribution of mass forces at the $k$-th sub-stage includes [30]:

equilibrium equations:

(5)$\mathrm{div}\hat{\mathsf{\sigma }}=0, x\in V_k,$

where $\hat{\mathsf{\sigma }}\left(x,t\right)$ is the stress tensor;

geometric Cauchy relations:

(6)$\hat{\mathsf{\epsilon }}=\frac{1}{2}\left(\nabla u+{\left(\nabla u\right)}^T\right), x\in V_k,$

where $u\left(x,t\right)$ is the displacement vector and $\hat{\mathsf{\epsilon }}\left(x,t\right)$ is the tensor of total deformations;

boundary conditions in displacements:

(7) $u=U, x\in S_{u,k},$

and tensions:

(8)$\hat{\mathsf{\sigma }}\cdot n=P, x\in S_{\mathsf{\sigma },k},$

where $S_u$ and $S_{\sigma }$ are the parts of the boundary with the specified displacements and loads, respectively.

Thermomechanical parameters of the material in the zone of “killed” elements exclude physical nonlinearity and are ideally elastic with degraded values:

${}_{}{}^4\widehat{C}\left(x\right), x\in V_k^{kil}<<{}_{}{}^4\widehat{C}\left(x,T\right), x\in V_k^{liv},$

where ${}_{}{}^4\widehat{C}$ is the fourth-rank tensor of elastic constants of the material.

2.2. Constitutive Relations for the PEEK Material and Their Adaptation to the Physical Models Available in ANSYS

The general system of equations of the boundary value problem of solid mechanics also includes constitutive relations.

The viscoelastic Prony model using the sum of exponentials with a constant volume compression modulus as the relaxation core is chosen as the constitutive relations:

(9)$\mathsf{\sigma }={{\displaystyle \int }}_0^{t}2G\left(t-\mathsf{\tau }\right)\frac{de}{d\tau }d\mathsf{\tau }+IK\theta ,$

where $\mathsf{\sigma }$ is the Cauchy stress tensor, $e$ is the deviatory part of the deformations, $I$ is the unit tensor, $K$ is the volume compression modulus, and $G\left(t\right)$ is the shear modulus:

(10)$G\left(t\right)=G_0\left[{\mathsf{\alpha }}_{\infty }^G+{\displaystyle \sum }_{i=1}^{n_G}{\mathsf{\alpha }}_i^G\exp \left(-\frac{t}{{\mathsf{\tau }}_i^G}\right)\right],$

where ${\mathsf{\alpha }}_i^G$ is the relative shear modulus for the shear relaxation times ${\mathsf{\tau }}_i^G$ and $n_G$ is the number of shear relaxation times. From the conditions $G_0={G|}_{t=0}$, $G_{\infty }={G|}_{t=\infty }$, it follows that:

(11)${\mathsf{\alpha }}_{\infty }^G=\frac{G_{\infty }}{G_0}, {\displaystyle \sum }_{i=1}^{n_G}{\mathsf{\alpha }}_i^G=\frac{G_0-G_{\infty }}{G_0}.$

2.2.1. Experimental Determination of Model Parameters

The use of the viscoelastic Prony model implies a preliminary determination of the material constants of the relaxation function. It is assumed that the material experiences only shear relaxation, and the volume compression modulus is constant in the model presented above. Traditionally, data from uniaxial tension-compression tests, rather than tests for shear, are used to obtain experimental parameters. The relaxation module for uniaxial stretching–compression has the form similar to (10):

$E\left(t\right)=E_0\left[c_{\infty }+{\displaystyle \sum }_{i=1}^{N_e}{c}_i\exp \left(-\frac{t}{{\mathsf{\beta }}_i}\right)\right],$

where $c_i$ are the relative modules of tension–compression for relaxation times ${\beta }_i$ and $N_e$ is the number of relaxation times of tension–compression. From the conditions $E_0={E|}_{t=0}$, $E_{\infty }={E|}_{t=\infty }$, it follows that:

(12)$c_{\infty }=\frac{E_{\infty }}{E_0}, {\displaystyle \sum }_{i=1}^{N_e}{c}_i=\frac{E_0-E_{\infty }}{E_0}.$

Assuming that the relaxation times for shear (${\tau }_i^G$) and stretch–compression (${\beta }_i$) coincide and their number is equal to ($N_e=n_G=n$). Then, from (11) and (12), we can obtain a relation of the form:

${\displaystyle \sum }_{i=1}^n{\mathsf{\alpha }}_i^G={\displaystyle \sum }_{i=1}^n{c}_i\left[\frac{G_0-G_{\infty }}{G_0}\cdot \frac{E_0}{E_0-E_{\infty }}\right].$

Denoting $B=\left[\frac{G_0-G_{\infty }}{G_0}\cdot \frac{E_0}{E_0-E_{\infty }}\right]$, we receive

(13)${\mathsf{\alpha }}_i^G=c_i\cdot B,$

It is known that

(14)$G_0=\frac{E_0}{2\left(1+{\nu }_0\right)},$

(15)$G_{\infty }=\frac{E_{\infty }}{2\left(1+{\nu }_{\infty }\right)},$

where ${\nu }_0$ and ${\nu }_{\infty }$ are the values of the Poisson’s ratio at the times $t=0$ and $t=\infty$, respectively. We find ${\nu }_{\infty }$ from the condition of constancy of the volume compression modulus $K^{*}$:

(16)$K^{*}=\frac{E_0}{3\left(1-2{\nu }_0\right)}=const,$

(17)$K_{\infty }=K^{*}=\frac{E_{\infty }}{3\left(1-2{\nu }_{\infty }\right)}\Rightarrow {\nu }_{\infty }=0.5\left(1-\frac{E_{\infty }}{3K^{*}}\right).$

Substituting (17) into (15), we receive:

(18)$G_{\infty }=\frac{E_{\infty }}{2\left(1+0.5\left(1-\frac{E_{\infty }}{3K^{*}}\right)\right)}=E_{\infty }\frac{3}{9-\frac{E_{\infty }}{K^{*}}}.$

Consequently, from (13), (16), (14), and (18) it is possible to calculate ${\mathsf{\alpha }}_i^G$, $K^{*}$, $G_0$, and $G_{\infty }$, having previously found $c_i$, $E_0$, $E_{\mathrm{\infty }}$, and ${\nu }_0$ from the tension–compression experiment. For this, we used the constitutive relations for the uniaxial case in the form:

(19)$\mathsf{\sigma }\left(t\right)={{\displaystyle \int }}_0^t\left[E_{\infty }+E_0{{\displaystyle \sum }}_{i=1}^{N_e}{c}_i\exp \left(-\frac{t-\mathsf{\tau }}{{\mathsf{\beta }}_i^{\prime }}\right)\right]d\mathsf{\epsilon }\left(\mathsf{\tau }\right),$

where ${\mathsf{\beta }}_i^{\prime }=\frac{{\mathsf{\beta }}_i}{A\left(T\right)}$ is the reduced time, $A\left(T\right)$ is the shift function. The material is assumed to be thermorheologically simple, so the Williams–Landel–Ferry shift function is used:

(20)$\lg \left(A\left(T\right)\right)=\frac{C_1\left(T-T_r\right)}{C_2+\left(T-T_r\right)},$

where $T$ is the current temperature, $T_r$ is the constant base temperature, and $C_1$ and $C_2$ are the empirical constants for the material. When representing the deformation as a harmonic function:

$\mathsf{\epsilon }\left(t\right)={\mathsf{\epsilon }}_a\sin \left(\omega t\right), (\omega =\mathit{const},T=\mathrm{var}), t\to \infty$

expression (19) is converted to the form:

$\mathsf{\sigma }\left(t\right)={\mathsf{\epsilon }}_a\{\sin \left(\omega t\right)[E_{\infty }+{\omega }^2{E}_0{\displaystyle \sum }_{i=1}^{n_G}\frac{{\mathsf{\beta }}_i^{\prime }{}^2{c}_i}{1+{\mathsf{\beta }}_i^{\prime }{}^2{\omega }^2}]+\cos \left(\omega t\right)\omega E_0{\displaystyle \sum }_{i=1}^{n_G}\frac{{\mathsf{\beta }}_i^{\prime }{c}_i}{1+{\mathsf{\beta }}_i^{\prime }{}^2{\omega }^2}\}.$

Hence, the expressions for the real and imaginary parts of the complex modulus will have the form:

(21)$E^{\prime }=E^{\infty }+{\omega }^2{E}_0{\displaystyle \sum }_{i=1}^{n_G}\frac{{\mathsf{\beta }}_i^{\prime }{}^2{c}_i}{1+{\mathsf{\beta }}_i^{\prime }{}^2{\omega }^2},$

(22)$E^{″}=\omega E_0{\displaystyle \sum }_{i=1}^{n_G}\frac{{\mathsf{\beta }}_i^{\prime }{c}_i}{1+{\mathsf{\beta }}_i^{\prime }{}^2{\omega }^2},$

Thus, Equations (9)–(20) constitute the PEEK physical model (Prony model and WLF), the parameters (material constants of the relaxation function and WLF constants) of which can be calculated from (21) and (22), knowing the temperature dependence of the real and imaginary parts of the PEEK complex modulus.

2.2.2. Identification of PEEK Properties according to Third-Party Sources Taking into Account the Kinetics of Crystallization

The numerical values of the temperature dependence function of the complex modulus firstly were determined from experimental data of third-party sources [7,16,17,18]. The papers describe (1) PEEK crystallization kinetics taking into account the cooling rate, (2) the PEEK viscoelastic model for predicting thermomechanical properties, and (3) the PEEK thermal deformation model.

Crystallization model formulation. The crystallization kinetics in semi-crystalline thermoplastic composites is considered in connection with the influence of the degree of crystallinity on the mechanical properties and possible contribution of volumetric shrinkage deformation. For a pure PEEK in dynamic conditions, it is set as follows:

(23)$X_{vc}=X^{\infty }{}_{vc}\left(w_1{F}_{vc1}+w_2{F}_{vc2}\right),$

with weight coefficients $w_1$ and $w_2=1-w_1$, where

(24)$$\begin{gathered} F_{vc1}=1-\exp \left[-C_{11}{{\displaystyle \int }}_0^{t}T\exp \left\{-n_1{t}^{\left(n_1-1\right)}\left[\frac{C_{21}}{T-T_g+51.6}+\frac{C_{31}}{T{\left(T_{m1}-T\right)}^2}\right]\right\}dt\right], \\ {F}_{vc2}=1-\exp \left[-C_{12}{{\displaystyle \int }}_0^{t}T\exp \left\{-n_2{t}^{\left(n_2-1\right)}\left[\frac{C_{22}}{T-T_g+51.6}+\frac{C_{32}}{T{\left(T_{m2}-T\right)}^2}\right]\right\}dt\right], \end{gathered}$$

where $X_{vc}$ is the volume fraction of the crystal phase, $X^{\infty }{}_{vc}$ is the equilibrium volume fraction, $T_g$ is the glass transition temperature, $T_{mi}$ is the melting point of crystals for the process ($i=1,2$), and $C_{1i}$, $C_{2i}$, and $C_{3i}$ are constants.

Viscoelastic model formulation. The real and imaginary parts of compliance of the material $J^{\prime }$ and $J^{″}$ are composed of amorphous ($J_{am}^{\prime }$, $J_{am}^{″}$ ) and crystalline ($J_{cr}^{\prime }$, $J_{cr}^{″}$) components:

(25)$$\begin{gathered} J^{\prime }=J_{am}^{\prime }\left(1-X_{vc}\right)+J_{cr}^{\prime }{X}_{vc}, \\ {J}^{″}=J_{am}^{″}\left(1-X_{vc}\right)+J_{cr}^{″}{X}_{vc}. \end{gathered}$$

The thermorelaxation transition is described by the following relations:

(26)$$\begin{gathered} J_{am}^{\prime }=J_{ua}+\left(J_{ua}-J_{ra}\right)+\left(\cos \mathsf{\psi }{)}^{\mathsf{\alpha }}\cos (\mathsf{\alpha }\mathsf{\psi }\right), \\ {J}_{am}^{″}=\left(J_{ra}-J_{ua}\right)\left(\cos \mathsf{\psi }{)}^{\mathsf{\alpha }}\sin (\mathsf{\alpha }\mathsf{\psi }\right). \end{gathered}$$

The real part of the compliance of the crystalline phase of the material is assumed to be constant, and the imaginary part is equal to zero:

(27)$$\begin{gathered} J_{cr}^{\prime }=J_{uc}, \\ {J}_{cr}^{″}=0, \end{gathered}$$

(28)$\mathsf{\psi }=\mathrm{arctg}\left({\omega \mathsf{\tau }}_{am}\right)$

where $J_{ua}$ is the compliance of the amorphous phase in the glassy state (instantaneous), $J_{ra}$ is the compliance of the amorphous phase in the equilibrium state (long-term), $J_{uc}$ is the elastic compliance of the crystalline phase, $\omega$ is the angular frequency of deformation of the sample in the DMA experiment, and $\mathsf{\alpha }$ is the relaxation distribution coefficient (in the range from 0 to 1), taking into account the tilt shift of the dynamic mechanical module in the glass transition domain due to changes in the degree of crystallinity.

The relaxation time of the amorphous phase ${\mathsf{\tau }}_{am}$ is directly proportional to the viscosity $\mathsf{\mu }$ in the model of a standard viscoelastic body:

(29)${\mathsf{\tau }}_{am}=\mathsf{\mu }\left(J_{ra}-J_{ua}\right).$

The temperature dependence of the delay time is described by two different laws: Arrhenius and Williams–Landel–Ferry (WLF). The first one operates in the range below the glass transition temperature, while:

(30)$\mathsf{\tau }={\mathsf{\tau }}_0{e}^{\frac{E_a\left(1/T-1/T_0\right)}{R}} \mathrm{for} T\le T_{0a}.$

The base temperature $T_0$ is assumed to be a linear function of the degree of crystallization in the range from $T_{0a}$ to $T_{0c}$. $E_a$ and $R$ are the specific activation energy and the universal gas constant, respectively. The approximation of WLF is given by the expression:

(31)$\mathsf{\tau }={\mathsf{\tau }}_0{10}^{\frac{-C_1\left(T-T_0\right)}{C_2+\left(T-T_0\right)}} \mathrm{for} T>T_{0a}.$

Constants $C_1$ and $C_2$ are determined from experiments on deformation at variable temperatures. The relaxation experiments are carried out in the interval up to 170 °C to fully identify the model (30, 31).

From the relations below, it is possible to calculate the real and imaginary parts of the complex modulus of the material:

(32)$$\begin{gathered} E^{\prime }=\frac{J^{\prime }}{{J^{\prime }}^2+{J^{″}}^2}, \\ {E}^{″}=\frac{J^{″}}{{J^{\prime }}^2+{J^{″}}^2}. \end{gathered}$$

Thermal deformation model formulation. A model based on the following hypotheses was chosen to describe the temperature deformation of the PEEK material: (1) the coefficient of thermal expansion of the material is constant over the temperature ranges of the glassy and highly elastic states. At the same time, it takes the following values:

(33) $$\begin{gathered} {\mathsf{\alpha }}^{th}\left(T\right)={\mathsf{\alpha }}_1,\hspace{0.17em}T<T_{g1}\hspace{0.17em}\mathit{glassy\; state} \\ {\mathsf{\alpha }}^{th}\left(T\right)={\mathsf{\alpha }}_2,\hspace{0.17em}T>T_{g2}\hspace{0.17em}\mathit{highly\; elastic\; state} \end{gathered}$$

The boundaries of the glass transition interval depend on the experimentally observed width of the relaxation transition:

$$\begin{gathered} T_{g1}=T_g-\mathsf{\Delta }{T}_g/2, \\ {T}_{g2}=T_g+\mathsf{\Delta }{T}_g/2, \end{gathered}$$

where $T_g$, $\mathsf{\Delta }{T}_g$ are the temperature and the glass transition interval, respectively. (2) In the interval of the relaxation transition, the coefficient of linear thermal expansion (CLTE) is calculated by the formula:

(34)${\mathsf{\alpha }}^{th}\left(T\right)=\frac{{\mathsf{\alpha }}_2-{\mathsf{\alpha }}_1}{2}\sin \left(\mathsf{\pi }\frac{T-T_g}{\mathsf{\Delta }{T}_g}\right)+\frac{{\mathsf{\alpha }}_2+{\mathsf{\alpha }}_1}{2},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{T}_{g1}\le T\le T_{g2}{.}_{}$

At the same time, the CTLE is interpreted as follows:

(35)${\mathsf{\alpha }}^{th}\left(T\right)=\frac{d{\mathsf{\epsilon }}^{th}\left(T\right)}{dT}$

(36)${\mathsf{\epsilon }}^{th}\left(T\right)={{\displaystyle \int }}_{T_0}^T{\mathsf{\alpha }}^{th}\left(T\right)dT$

Total free deformation also includes the shrinkage of the material ${\mathsf{\epsilon }}^{cr}$ as a result of the partial crystallization of PEEK:

(37)${\mathsf{\epsilon }}^f={\mathsf{\epsilon }}^{cr}+{\mathsf{\epsilon }}^{th}.$

The expression for ${\mathsf{\epsilon }}^{cr}$ can be obtained by knowing the dependence of the density of the material on the volume content of the crystalline phase:

$\rho \left(X_{vc}\right)={\rho }_{am}\left(1-X_{vc}\right)+{\rho }_{cr}{X}_{vc},$

herewith

(38)${\mathsf{\epsilon }}^{cr}\left(X_{vc}\right)=1/3\frac{\left({\rho }_{am}-\rho \left(X_{vc}\right)\right)}{{\rho }_{am}}.$

Thus, Equations (23)–(32) constitute the PEEK physical model (micromechanical model of crystallization kinetics taking into account the cooling rate), and Equations (33)–(38), the PEEK thermal deformation model. These expressions were used to reproduce the experimentally obtained graphs [7], which, in turn, were used to identify the constants of the Prony and WLF models. It was carried out as follows.

Crystallization model calculation. The identification of the parameters of the crystallization model (23-24) was carried out in [16]. The verification showed that the parameter values presented in the paper (found from the condition of the minimum discrepancy between the calculation and the experiment) fail to ensure the achievement of the claimed accuracy in practice. Therefore, the search for their optimal values was implemented in the Matlab package. The following vector of the required constants varied:

$x=\left\{X_{vc\infty },n_1,n_2,T_{m1},T_{m2},C_{11},C_{12},C_{21},C_{22},C_{31},C_{32},T_g,w_1\left({\dot{T}}_i\right)\right\},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$

The method of a deformable polyhedron (Nelder–Mead) was used to search for the minimum specific discrepancy of the calculated and experimental data:

(39)$\mathsf{\Phi }\left(x\right)=\frac{{{\displaystyle \sum }}_{i=1}^{M_e}\left|X_{vci}^e-X_{vci}^s\left(x\right)\right|}{{{\displaystyle \sum }}_{i=1}^{M_e}\left|X_{vci}^e\right|}\to min.$

The experimental data of paper [7], presented in Figure 2 by red curves, were used as reference data. The experimental curves of the degree of crystallization for cooling rates $T_i$ = 1, 10 and 35 K·s⁻¹ are given in the work.

The result of (39) is presented in Table 1 and Table 2.

Viscoelastic model calculation. The results from calculating the temperature dependence of the complex modulus according to (32) from paper [7] were used as a reference to identify the parameters of (9)–(20) available in ANSYS. The values of the material constants accepted in the calculation are presented in Table 3.

The calculation result is shown in Figure 3a,b shows experimental graphs of the complex modulus amplitude from [7].

Thermal deformation model calculation. Figure 4 shows the results of calculating the free deformation ${\mathsf{\epsilon }}^f={\mathsf{\epsilon }}^{cr}+{\mathsf{\epsilon }}^{th}$ and its component obtained from (33)–(38). The thermomechanical constant values of the model are presented in Table 4.

Identification of the Prony model constants from graph data. The temperature dependence of the complex modulus at the cooling rate of $\dot{T}$= 1 K·s⁻¹, shown in Figure 3a with thin lines, was used to determine the parameters of the Prony model with the temperature–time analogy WLF (9)–(20). The formulation of the problem of finding constants of the model is similar to (39). The obtained values are shown in Figure 5 and Figure 6. Figure 5 shows the temperature dependence of the real $E^{\prime }$ and imaginary $E^{″}$ parts of the complex modulus obtained during the experiment and numerical solution. Table 5 and Table 6 show the instantaneous properties of the material and the temperature–time shift constants, respectively. The dependence of the weight coefficients $c_i$ on the relaxation times ${\mathsf{\beta }}_i$ is shown in Table 7 and Figure 6.

Figure 7 shows the result of testing the model (9)–(20) in the finite element package of ANSYS. The calculation was performed without taking into account temperature and shrinkage deformations. At the initial moment, the virtual cubic sample instantly (in 10⁻⁹ s) is stretched along the x-axis by an amount of ${\mathsf{\epsilon }}_x$ = 0.01. At the same time, the faces normal to the other axes remain free. Thus, a uniaxial stress state (USS) is realized. Further, the sample is heated at a fixed deformation at a constant rate from room temperature to 200 °C. The transition to a highly elastic state is accompanied by a drop in modulus and voltage. It can be seen from the figure that the transition zone shifts to a zone of higher temperatures with an increase in the heating rate, which is typical for viscoelastic polymers of the studied class. Thus, the adequacy of the obtained model of the thermomechanical behavior of the PEEK material is qualitatively confirmed.

2.2.3. Identification of PEEK Properties Based on the Results of Our DMA Experiment

When using the data obtained from [7,16,17,18] in numerical simulation, it was not possible to obtain results close to the full-scale experiment. The resulting deflections of PEEK plates differed significantly from the experimental ones even with an increase in the degree of detail of the slicer’s trajectory. Due to the fact that the data presented in these papers were obtained at room temperature under loading conditions that are not typical for the printing process, it was decided to conduct our own DMA experiment in a wide temperature range.

To determine the numerical values of the temperature dependence function of the complex modulus, full-scale and numerical experiments were carried out, and the model was verified. A thermomechanical experiment for the glass-transiting material was carried out using a dynamic mechanical analyzer DMA Q800 TA.

After separation from the platform, the samples showed residual bending, the direction of which (edges up from the platform) can be seen in Figure 8. Next, the samples were subjected to heat treatment according to two different scenarios:

Scenario 1. Heating in the furnace up to 250 °C, followed by cooling together with the furnace (approximate cooling rate 2760 °C/min = 46 °C/s);

Scenario 2. Heating in the furnace to 250 °C, followed by cooling in running water (approximate cooling rate 3.8 °C/min = 0.064 °C/s). Only 2 samples.

It was assumed that for samples cooled at such different rates, it would be possible to reveal differences in thermomechanical behavior due to different conditions for the formation of the crystal structure of the material and the final values of the degree of crystallinity (see Section 3). Only 2 samples.

It was observed that the residual bending after additional heat treatment for all samples increased by approximately 1.5 times.

A three-point bending flexural test was carried out in accordance with the requirements of State All-Union standard 4648-2014 [31], according to which, the test samples were formed with overall dimensions: 4.7 × 1.76 × 79 mm (Figure 1).

Next, 2 experiments were carried out on a three-point bending for a couple of “fast”–“slow” samples:

1. In the temperature range 30–220 °C;

2. In the temperature range 30–330 °C. Figure 8 shows a sample after these tests.

In both cases, the frequency of the kinematic impact was 1 Hz.

The results of measuring the parameters of the complex modulus are shown in Figure 9. From the figure, in particular, it follows that there is no effect of the cooling rate of the sample on its thermomechanical behavior. Graphs marked with a dash correspond to cooling in water, and solid graphs correspond to cooling in air. The stiffness of the “fast” sample in Figure 9a is less, and in Figure 9b, it is greater than that of the one cooled together with the furnace. In addition, there are no significant shifts of the drop region of the complex modulus $E^{\prime }$ for different samples up or down the temperature scale. We can conclude that the differences in the experimental curves lie within the limits of the experimental error.

Identification of the Prony model constants from graph data. The search for unknown constants was carried out by the method of nonlinear programming in the MatLab package using the fminsearch function. The problem of finding the minimum of the combined discrepancy between the experimental characteristics and those calculated by the model (21) and (22) was solved:

(40)$\mathsf{\Phi }\left(x\right)={\mathsf{\Phi }}_{E^{\prime }}\left(x\right){\mathsf{\Phi }}_{E^{″}}\left(x\right)\to \min ,$

where ${\mathsf{\Phi }}_{E^{\prime }}\left(x\right)={\displaystyle \sum }_{j=1}^{N_g}\left|{E^{\prime }}_j^{\exp }-{E^{\prime }}_j^{sol}\right|/N_g$, ${\mathsf{\Phi }}_{E^{″}}\left(x\right)={\displaystyle \sum }_{j=1}^{N_g}\left|{E^{″}}_j^{\exp }-{E^{″}}_j^{\mathit{sol}}\right|/N_g$, and $sol$ are the experimental and computational values at each point of measurement. Additional limitations were put on the non-negativeness of the coefficients ${\mathsf{\alpha }}_i^G$: ${\mathsf{\alpha }}_i^G\ge 0$.

The data obtained after solving Equation (40) are given in Table 8, Table 9 and Table 10.

Table 8

Instantaneous properties of viscoelastic material.

${\mathit{E}}_0$, Pa	${\mathit{\nu }}_0$
$3.0145420\times {10}^9$	0.42

Table 9

Values of constants of the Williams–Landel–Ferry shift function.

${\mathit{T}}_{\mathit{r}}$ , K	${\mathit{C}}_1$	${\mathit{C}}_2$
318	84.01	542.8

Table 10

Values of viscoelastic constants of the relaxation function for 30 relaxation times.

$\mathit{i}$	1	2	3	4	5	6	7	8
${\mathsf{\alpha }}_{\mathit{i}}^{\mathit{G}}$	$9.782\times {10}^{-3}$	$1.111\times {10}^{-2}$	$1.240\times {10}^{-2}$	$7.130\times {10}^{-3}$	$9.511\times {10}^{-3}$	$1.061\times {10}^{-2}$	$1.121\times {10}^{-2}$	$1.193\times {10}^{-2}$
${\beta }_{\mathit{i}}$	$1.636\times {10}^{-1}$	$2.575$	$4.053\times 10$	$6.379\times {10}^2$	$1.004\times {10}^4$	$1.458\times {10}^5$	$2.488\times {10}^6$	$3.915\times {10}^7$
$\mathit{i}$	9	10	11	12	13	14	15	16
${\mathsf{\alpha }}_{\mathit{i}}^{\mathit{G}}$	$1.340\times {10}^{-2}$	$1.556\times {10}^{-2}$	$1.893\times {10}^{-2}$	$7.206\times {10}^{-2}$	$3.082\times {10}^{-1}$	$2.394\times {10}^{-1}$	$9.227\times {10}^{-2}$	$3.198\times {10}^{-2}$
${\beta }_{\mathit{i}}$	$6.163\times {10}^8$	$9.701\times {10}^9$	$1.527\times {10}^{11}$	$2.403\times {10}^{12}$	$3.783\times {10}^{13}$	$5.954\times {10}^{14}$	$9.372\times {10}^{15}$	$1.475\times {10}^{17}$
$\mathit{i}$	17	18	19	20	21	22	23	24
${\mathsf{\alpha }}_{\mathit{i}}^{\mathit{G}}$	$1.551\times {10}^{-2}$	$1.088\times {10}^{-2}$	$9.063\times {10}^{-3}$	$8.216\times {10}^{-3}$	$6.637\times {10}^{-3}$	$5.990\times {10}^{-3}$	$6.824\times {10}^{-3}$	$8.466\times {10}^{-3}$
${\mathsf{\beta }}_{\mathit{i}}$	$2.322\times {10}^{18}$	$3.655\times {10}^{19}$	$5.753\times {10}^{20}$	$9.055\times {10}^{21}$	$1.425\times {10}^{23}$	$2.243\times {10}^{24}$	$3.531\times {10}^{25}$	$5.558\times {10}^{26}$
$\mathit{i}$	25	26	27	28	29	30
${\mathsf{\alpha }}_{\mathit{i}}^{\mathit{G}}$	$8.661\times {10}^{-3}$	$5.874\times {10}^{-3}$	$4.392\times {10}^{-3}$	$5.120\times {10}^{-3}$	$7.408\times {10}^{-3}$	$8.428\times {10}^{-3}$
${\beta }_{\mathbf{i}}$	$8.749\times {10}^{27}$	$1.377\times {10}^{29}$	$2.167\times {10}^{30}$	$3.412\times {10}^{31}$	$5.370\times {10}^{32}$	$8.452\times {10}^{33}$

Distribution of weight coefficients $c_i$ according to corresponding relaxation times ${\mathsf{\beta }}_i$ is presented in Figure 10.

Figure 11a,b shows a comparison of the obtained solution of the minimization problem (40) with the experiment. It should be noted that the obtained model is suitable for describing the thermomechanical behavior.

2.3. Numerical Algorithm

2.3.1. Solution of Thermal Conductivity Problem

The algorithm used to calculate temperature fields in the ANSYS finite element package in the numerical simulation of the layer-by-layer deposition involves the following calculation steps:

1. The creation of a finite element model, including the volume of the future product and the platform with the corresponding thermophysical properties (Figure 12a). The thermophysical properties of the material are taken from paper [7]. At the same time, in accordance with the data presented in the paper, thermal conductivity $\mathsf{\lambda }$ = 0.251 W/(m K) and heat capacity $C_p$ = 1399 J/(kg K) were assumed to be constant. For the platform, it is Steel;

2. “Killing” (EKILL command) the elements of the product that is absent in the actual deposition process before it starts (Figure 12b);

3. The heating of the platform up to 140 °C (413.15 K) according to Formula (3) and waiting for 300 s for the thermal state of the platform to reach the stationary mode (Figure 12c);

4. In the cycle by the sub-stages of the passage of the deposition zones of the next calculated section of the plate material:

4.1. “Aliving” (EALIVE command) of all elements of the k-th deposition zone (Figure 12d);

4.2. Setting the conditions of convective heat transfer at the free surfaces of the platform and “aliving” the plate material according to Formula (2) (SF command) (Figure 12e);

4.3. The instantaneous heating up to 430 °C (703.15 K) of the k-th deposition zone elements by a source of thermal energy distributed over the volume (see Formula (1)) (Figure 12f). The appropriate amount of energy that the extruder heater takes for heating PEEK up to 430 °C is supplied to the model by instantaneously applying the appropriate temperature throughout the volume of the k-th deposition zone of the PEEK plate material using the D command;

4.4. Removing the heat source and waiting for the partial cooling of the system for a time interval (waiting stage) equal to the waiting time $t_w$ (Figure 12g).

The waiting time $t_w$ depends on the degree of detail of the trajectory of the heat source and is calculated by Formula (41):

(41)$t_w=\frac{l}{v\cdot n_{regX}\cdot n_{regZ}},$

where $l$ is the plate length, $n_{regX}$ and $n_{regZ}$ are the number of sections of the plate material into which the trajectory of the heat source is divided along the long and short sides of the sample, respectively, and $v$ is the printing speed.

2.3.2. Solving the Problem of Determining the Stress–Strain State of a Structure

The algorithm is similar to that in Section 2.3.1.

1. The creation of a finite element model, including the volume of the future product and the platform with the corresponding thermomechanical properties (Figure 13a). The following values of PEEK parameters were used in ANSYS. The values of the constants from Table 8, Table 9 and Table 10 were used for its viscoelastic model, and the data from Figure 4 were used for the simulation of its temperature deformations:

1. “Killing” (EKILL command) of the elements of the product that are absent in the deposition process before it starts. Setting boundary conditions in displacements (D command) (Figure 13b);

2. In the cycle by sub-stages of the passage of the deposition zones of the next calculated section of the plate material:

2.1. “Alive” (EALIVE command) all elements of the k-th deposition zone (Figure 13c);

2.2. The application of previously calculated temperatures for a given point in time to the h nodes of the model during the time of exposure to a heated filament (instantaneously) (LDREAD command) (Figure 13d);

2.3. The sequential application of previously calculated temperatures for the waiting stage during $t_w$ seconds (LDREAD command) (Figure 13e);

3. Ambient temperature application to the nodes of the model, VAT calculation (Figure 13f).

4. The separation of the plate from the platform and the calculation of the VAT (Figure 13g). After separation, the value of the residual bending was taken from the middle line. The data were obtained through the PATH command.

3. Results and Discussion

3.1. The Results of the Experimental Study of the Residual Bending of PEEK Samples Created Using Different Deposition Trajectories

The deposition of the samples was carried out using a standard FDM 3D printer equipped with a sealed thermal chamber and a high-temperature extruder (up to 500 °C) of our own production, the design and principle of operation of which are described in [32]. The material used was a wire/filament composed of PEEK with a diameter of 1.75 mm, the parameters of which are presented in [33], manufactured by the company CreatBot.

The sample was a rectangular parallelepiped with overall dimensions: 4.7 × 1.76 × 79 mm. The samples were obtained by the method of additively manufactured (FDM) with a longitudinal and transverse spiral layout of beads 0.71 mm wide and 0.29 mm high (Figure 14). In total, 5 samples were created with a longitudinal layout of the material and 5 with a transverse one. The printing speed was 15 mm/s, the ambient temperature was 100 °C, and the platform heating temperature was 140 °C. The test samples were printed with a software-specified 100% filling of the internal volume.

The results were processed as follows: (1) the samples were placed in a plane perpendicular to the optical axis of the system and photographed; (2) the bending was digitized in the GetData Graph Digitizer program by specifying points on the image (Figure 15a); (3) according to the obtained data, bending graphs were plotted in Matlab (Figure 15b), which were then aligned relative to the extreme points (Figure 15c), and the converted curves to symmetrical (Figure 15d); (4) values at coinciding points were averaged and approximated (Figure 15e).

3.2. Solving the Problem of Finding Stress–Strain State of PEEK Products

The values of the material constants obtained from works [7,16,17,18] when they were used in the numerical simulation did not allow obtaining results close to the experiment. The numerical deflection was very different from the results of a full-scale experiment, regardless of the degree of detail of the trajectory of the heat source and the improvement of the numerical algorithm (Section 3.3).

As it was mentioned before, we assumed that this problem is due to the fact that the data presented in these papers were obtained at room temperature under loading conditions that are not typical for the printing process. As a result, it was decided to conduct our own experiment on DMA in a wide temperature range (Section 2.2.3). The use of the material constants of the relaxation function obtained in our own experiment helped to obtain acceptable results. A numerical simulation with new constants obtained from our own experiment was carried out on the same model.

Figure 16a,b shows the bendings obtained in three numerical solutions (which differ in the degree of discretization of the trajectory of the heat source (Section 3.3)) and a full-scale experiment. The results show a good agreement between the experimental data and the data obtained from the simulation.

A significant error is observed only in the roughest version of the calculation (when the layers with an area equal to the longitudinal cross-section of the sample were deposited). When solving the problem of thermomechanics, only those temperature fields are read that correspond to the moments of deposition of a new section of the plate material and its cooling by the time of the deposition of the next one; the values between them are linearly interpolated. There are only 12 such solution sub-steps in a rough calculation, which, as we assumed, may not be enough for the correct application of temperatures. However, an increase in the number of sub-steps only led to a slight increase in the accuracy of the solution, from which it can be concluded that the roughest version has a great error due to the fact that this method of applying the material is not suitable for modeling an additive process. Further, the roughest variant is not considered.

3.3. Numerical Solution Features

Adhesion of the sample to the platform during the deposition and its separation upon completion of the structure formation was modeled in two ways: (1) by linking joint knots in the contact zone (CPINTF command); (2) by applying surface-to-surface contact elements, TARGE170 and CONTA173. The second approach showed the best accuracy of the solution, in particular, adequate heat distribution during the instantaneous application of the first layer of the product (Figure 17), which made it possible to take into account the adhesion force, the excess of which leads to the premature detachment of the sample from the platform.

Additionally, the influence of the degree of discretization of the trajectory of the heat source generated by the slicer program on the quality of the solution was studied. The tasks were solved three times: in the first case, the samples were formed in layers with the area of the longitudinal cross-section of the sample; in the second case, these layers were applied with longitudinal beads; in the third case, each roller was assembled from even smaller sections of the plate material (Figure 18). To assess the proportionality and correctness of setting the calculated time and boundary conditions in three solutions, the energy balance was analyzed. The amount of heat of the system was calculated for the moment of exposure to the heating of the next calculated section of the plate material and the moment of the end of the waiting time according to the formula:

(42)$Q={\displaystyle \sum }_{i=1}^{n_e}{c}_i{\rho }_i{V}_i{T}_i,$

where $n_e$ is the number of “alive” elements of the system at the $k$-th sub-stage of the solution, and $V_i$, $c_i$, ${\rho }_i$, and $T_i$ are the volume, heat capacity, density, and average temperature of the $i$-th element. The energy balance was assessed sequentially. At first, the process was modeled without taking into account the heat exchange of the sample with the platform and the surrounding medium, which made it possible to correct the deposition time for sections of the plate material (Figure 19). When the heat exchange between the sample and the platform is included in the calculation (there is no heat exchange with the surrounding medium at this stage), excess energy began to accumulate in the system (Figure 20a). This was due to the fact that during the instantaneous heating of the $k$-th “alived” zone of the elements, heat was also transferred to the adjacent elements, which had common nodes with the “alive” zone. When using the 20-node elements for heating, only the middle nodes did not give the desired result: the “alived” section was heated unevenly, which was not true. The problem was solved by creating intermediate nodes near the common ones. Thus, when heating only the auxiliary nodes, the heat was evenly distributed over the volume of the “alived” section, did not spread to the adjacent elements (Figure 21), and was supplied in the required amount. When the remaining boundary conditions are taken into account, the energy balance was kept.

4. Conclusions

During the work, we established the fact that the data on the thermal and viscoelastic properties of PEEK presented in the literature are not quite adequate for additive manufacturing due to results presented being obtained at room temperature under loading conditions that are not typical for the deposition process. For this reason, our own DMA experiment was carried out to determine the numerical values of the temperature dependence function of the complex relaxation modulus in a wide temperature range. Using the ANSYS APDL language, our own parameterized algorithm was developed for calculating non-stationary temperature fields and the stress–strain state of PEEK products during their manufacturing using a wire-based deposition method (FDM). The influence of the degree of detailing of the heat source trajectory and the method of connecting the sample to the platform on the quality of the solution is analyzed. The heat supply mechanism was debugged, and the estimated time was coordinated with an increase in the detailing of the heat source trajectory by taking into account the heat balance. As part of the analysis of the results of the numerical simulation, a good agreement between the calculated data and the experiment was shown.

As a result of the study of the thermomechanics of the additive manufacturing of PEEK products, the following has been found:

1. The Prony viscoelastic model is suitable as a model for describing PEEK during hardening;

2. The literature data do not make it possible to numerically obtain results close to the full-scale experimental ones. The constants obtained from our own DMA experiment give a qualitative solution. In addition, the results of the DMA experiment demonstrate a smoother transition of the material from the highly elastic to the glassy state and indicate that the sample cooling rate does not affect its thermomechanical behavior;

3. Contact elements in the sample-platform zone minimize the calculated errors and allow taking into account the adhesion force in the contact zone; therefore, they are better suited for modeling sticking with the subsequent detachment of the sample from the platform;

4. The application of temperature to the common nodes of the adjacent sections of the plate material leads to the supply of excess heat to the system with the instantaneous heating of the “alived” section. You can fix this by creating additional nodes near the common ones.

5. Significant reduction in the detailing of the heat source trajectory leads to a great error in the solution.

At present, it is planned to carry out identification tests for shear and verification experiments on shell growth, which will allow us to assess the acceptability of the hypothesis of a constant bulk modulus.

Author Contributions

Conceptualization, D.N.T. and O.Y.S.; methodology, O.Y.S., A.A.A. and D.N.T.; software, O.Y.S. and A.A.A.; validation, O.Y.S.; formal analysis, O.Y.S. and A.A.A.; investigation, O.Y.S., A.A.A., A.A.O. and A.A.L.; resources, O.Y.S. and D.N.T.; data curation, O.Y.S., A.A.A. and A.A.L.; writing—original draft preparation, O.Y.S., A.A.A. and A.A.O.; writing—review and editing, A.A.O. and A.A.A.; visualization, O.Y.S., A.A.A. and A.A.L.; supervision, D.N.T. and O.Y.S.; project administration, D.N.T.; funding acquisition, D.N.T. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the fact that further study will be carried out using the same data.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Modeled geometry: the calculated scheme of the problem: (1) the domain occupied by the final product; (2) the platform on which deposition is carried out; [Forumla omitted. See PDF.]—the domain occupied by platform and part of the printed sample at the [Forumla omitted. See PDF.]-th sub-stage of the solution; [Forumla omitted. See PDF.]—free surfaces of the platform and part of the printed sample at the [Forumla omitted. See PDF.]-th sub-stage of the solution; [Forumla omitted. See PDF.]—heated platform surface.

Figure 2. Temperature dependence of the degree of crystallization of PEEK during cooling at different speeds: (a) without detailing the mechanisms of crystal nucleation; (b) with detailing. Black, green, and blue lines show calculations; red ones show experiments [7].

Figure 3. Temperature dependence of the complex PEEK module during cooling at different speeds: (a) calculation according to (32), where the blue lines are [Forumla omitted. See PDF.], red lines are [Forumla omitted. See PDF.], and black lines are amplitude [Forumla omitted. See PDF.]; (b) experimental data of the paper [7].

Figure 4. Temperature dependence of free deformation components of PEEK during cooling at different speeds. Black curves—[Forumla omitted. See PDF.], red—[Forumla omitted. See PDF.], blue—[Forumla omitted. See PDF.].

Figure 5. The result of the approximation of the temperature dependence of the complex PEEK module when cooled at a rate of [Forumla omitted. See PDF.] = 1 K·s−1: (a) is the real part; (b) is the imaginary part. “Experiment” is a reproduction of experiment results [7] by calculating (32), “Calculation” is the calculation by (21–22).

Figure 5. The result of the approximation of the temperature dependence of the complex PEEK module when cooled at a rate of [Forumla omitted. See PDF.] = 1 K·s−1: (a) is the real part; (b) is the imaginary part. “Experiment” is a reproduction of experiment results [7] by calculating (32), “Calculation” is the calculation by (21–22).

Figure 6. Dependence of [Forumla omitted. See PDF.] on [Forumla omitted. See PDF.].

Figure 7. Curves of the temperature–time relaxation of a pinched PEEK sample without taking into account free deformation (calculation in ANSYS).

Figure 8. PEEK sample before and after DMA experiment.

Figure 9. Results of DMA tests of PEEK samples cooled in water (dashed lines) and together with the furnace (solid lines).

Figure 9. Results of DMA tests of PEEK samples cooled in water (dashed lines) and together with the furnace (solid lines).

Figure 10. Distribution of weight coefficients [Forumla omitted. See PDF.] according to corresponding relaxation times [Forumla omitted. See PDF.].

Figure 11. Temperature dependence: (a) the real part of the complex modulus [Forumla omitted. See PDF.]. The thin line is the experiment; the thick line is the calculation; (b) the imaginary part of complex modulus [Forumla omitted. See PDF.]. The thin line shows the experiment, and the thick line shows the calculation.

Figure 11. Temperature dependence: (a) the real part of the complex modulus [Forumla omitted. See PDF.]. The thin line is the experiment; the thick line is the calculation; (b) the imaginary part of complex modulus [Forumla omitted. See PDF.]. The thin line shows the experiment, and the thick line shows the calculation.

Figure 12. Modeling in ANSYS by stages of thermal conductivity problem solution: (a) the finite element model; (b) calculated scheme after “killing” the elements; (c) initial conditions applied to the platform before the deposition of the plate material; (d) “aliving” of new zone of plate elements at the k-th stage of deposition (gray zone); (e) setting the conditions of convective heat transfer to system at the k-th stage of deposition; (f) instantaneous heating of the k-th deposition zone of plate elements; (g) removing the heat source and waiting for the partial cooling of the system.

Figure 13. Modeling in ANSYS by stages of thermomechanical problem solution: (a) the finite element model; (b) boundary conditions applied to the platform and “dead” elements of the plate, which are removed as soon as the k-th zone of the plate elements is “alive”; (c) “aliving” of the new zone of the plate elements at the k-th stage of deposition; (d) application of previously calculated temperatures after the heating of the zone of plate elements at the k-th stage of deposition; (e) application of previously calculated temperatures after the cooling of the zone of the plate elements at the k-th stage of deposition; (f) application of ambient temperature to system; (g) separation of the plate from the platform.

Figure 13. Modeling in ANSYS by stages of thermomechanical problem solution: (a) the finite element model; (b) boundary conditions applied to the platform and “dead” elements of the plate, which are removed as soon as the k-th zone of the plate elements is “alive”; (c) “aliving” of the new zone of the plate elements at the k-th stage of deposition; (d) application of previously calculated temperatures after the heating of the zone of plate elements at the k-th stage of deposition; (e) application of previously calculated temperatures after the cooling of the zone of the plate elements at the k-th stage of deposition; (f) application of ambient temperature to system; (g) separation of the plate from the platform.

Figure 14. General view (side and top) of the resulting product with the layout: (a) along; (b) across the long side.

Figure 15. Processing the experimental results: (a) digitizing the bending; (b) plotting graphs based on the data obtained; (c) aligning the curves with respect to the extreme points; (d) converting curves to symmetrical; (e) averaging the values at coinciding points with subsequent approximation. The graphs show the value of the deflection (w) along the length (l) of the sample.

Figure 15. Processing the experimental results: (a) digitizing the bending; (b) plotting graphs based on the data obtained; (c) aligning the curves with respect to the extreme points; (d) converting curves to symmetrical; (e) averaging the values at coinciding points with subsequent approximation. The graphs show the value of the deflection (w) along the length (l) of the sample.

Figure 16. Residual bending of the sample after its separation from the platform in three numerical solutions and a full-scale experiment: (a) for samples with a longitudinal laying of the material; (b) for samples with a transverse laying of the material. The graphs show the value of the deflection (w) along the length (l) of the sample.

Figure 16. Residual bending of the sample after its separation from the platform in three numerical solutions and a full-scale experiment: (a) for samples with a longitudinal laying of the material; (b) for samples with a transverse laying of the material. The graphs show the value of the deflection (w) along the length (l) of the sample.

Figure 17. System at [Forumla omitted. See PDF.] = [Forumla omitted. See PDF.] s: (a) application of command CPINTF; (b) application of contact elements TARGE170 and CONTA173.

Figure 18. Increasing the degree of detail of the heat source trajectory: (a) solid (1 step—1 layer); (b) beads (1 step—1/5 of the layer); (c) sections (1 step—1/40 of the layer).

Figure 19. Energy balance for solutions with varying degrees of discretization of the problem without taking into account the heat transfer of the sample with the platform and the surrounding medium: (a) throughout the entire deposition process; (b) during the formation of the first layer. Three hundred seconds is the waiting time for the thermal state of the platform to reach the stationary mode.

Figure 19. Energy balance for solutions with varying degrees of discretization of the problem without taking into account the heat transfer of the sample with the platform and the surrounding medium: (a) throughout the entire deposition process; (b) during the formation of the first layer. Three hundred seconds is the waiting time for the thermal state of the platform to reach the stationary mode.

Figure 20. Supply of energy to the system without heat exchange between the sample and the medium: (a) before correction; (b) after it. Three hundred seconds is the waiting time for the thermal state of the platform to reach the stationary mode.

Figure 21. Demonstration of the refiring of the mechanism of energy supply to the system: (a,b)—the time of the end of the cooling of the first layer and the instantaneous heating of the second layer of plate material, respectively, before correction; (c,d)—the same, after correction.

Table 1

Values of parameters of the PEEK material crystallization kinetics model.

${\mathit{X}}_{\mathit{v}\mathit{c}\mathit{\infty }}$	${\mathit{n}}_1$	${\mathit{n}}_2$	${\mathit{T}}_{\mathit{m}1}$ , K	${\mathit{T}}_{\mathit{m}2}$ , K	${\mathit{C}}_{11}$ , c−nK−1
0.2721	2.7986	1.4963	$5.9246\times {10}^2$	$6.1696\times {10}^2$	$3.4052\times {10}^{10}$
${\mathit{C}}_{12}$ , c−nK−1	${\mathit{C}}_{21}$ , K	${\mathit{C}}_{22}$ , K	${\mathit{C}}_{31}$ , K3	${\mathit{C}}_{32}$ , K3	${\mathit{T}}_{\mathit{g}}$ , K
$3.4052\times {10}^{10}$	$5.0698\times {10}^3$	$4.2589\times {10}^3$	$3.6516\times {10}^7$	$4.6533\times {10}^7$	$4.1700\times {10}^2$

Table 2

Weight coefficients of the mechanisms of crystallization of the PEEK material.

${\dot{\mathit{T}}}_{\mathit{i}}$ , K·s−1	1	10	35
${\mathit{w}}_{\mathbf{1}}\left({\dot{T}}_i\right)$	$8.1209\times {10}^{-2}$	$7.1103\times {10}^{-1}$	$9.9980\times {10}^{-1}$

Table 3

Values of constants for calculating the dependence of a complex modulus.

${\mathit{J}}_{\mathit{u}\mathit{a}}$, Pa−1	${\mathit{J}}_{\mathit{r}\mathit{a}}$, Pa−1	${\mathit{J}}_{\mathit{u}\mathit{c}}$, Pa−1	ω, rad s−1	τ0, s	Ea, J mol−1	R, J mol−1 K−1
$380\times {10}^{-12}$	$15000\times {10}^{-12}$	$120\times {10}^{-12}$	1	1	$600\times {10}^3$	8.314462618
T 0a , K	T 0c , K	${\mathit{C}}_1$	${\mathit{C}}_2$	$\alpha$	${\mathsf{\alpha }}_{\mathit{a}\mathit{m}}$	${\mathsf{\alpha }}_{\mathit{c}\mathit{r}}$
414	533	25	100	0.5	1	0.4

Table 4

Thermomechanical values of constants.

${\mathit{\rho }}_{\mathit{a}\mathit{m}}$, kg·m−3	${\mathit{\rho }}_{\mathit{c}\mathit{r}}$, kg·m−3	${\mathsf{\alpha }}_2$, K−1	${\mathsf{\alpha }}_2$, K−1	$\mathsf{\Delta }{\mathit{T}}_{\mathit{g}}$, K
$1262.6$	$1400.6$	$50\times {10}^{-6}$	$200\times {10}^{-6}$	10

Table 5

Values of instantaneous properties of viscoelastic material.

${\mathit{E}}_0$, Pa	${\mathit{\nu }}_0$
$3.2335877\times {10}^9$	0.3

Table 6

Values of constants of the Williams–Landel–Ferry shift function.

Tr, K	C 1	C 2
397.24846	61.682023	150.11259

Table 7

Values of viscoelastic constants of the relaxation function for 10 relaxation times.

i	1	2	3	4	5	6	7	8	9	10
${\mathsf{\alpha }}_{\mathit{i}}^{\mathit{G}}$	$1.37\times {10}^{-4}$	$2.62\times {10}^{-8}$	$8.92\times {10}^{-4}$	$1.10\times {10}^{-3}$	$3.78\times {10}^{-4}$	$6.29\times {10}^{-4}$	$1.12\times {10}^{-3}$	$1.05\times {10}^{-3}$	$4.96\times {10}^{-4}$	$9.74\times {10}^{-1}$
${\beta }_{\mathit{i}}$	$1.00\times {10}^{-4}$	$7.74\times {10}^{-3}$	$5.99\times {10}^{-1}$	$4.64\times 10$	$3.59\times {10}^3$	$2.78\times {10}^5$	$2.15\times {10}^7$	$1.67\times {10}^9$	$1.29\times {10}^{11}$	$1.00\times {10}^{13}$

References

1. Salamov, A.K.; Mikitaev, A.K.; Beev, A.A.; Beeva, D.A.; Kumysheva, Y.A. Polyetheretherketones (PEEK) as representatives of aromatic polyaryleneetherketones. Fundam. Res.; 2016; 16, pp. 63-66.

2. Ma, R.; Tang, T. Current Strategies to Improve the Bioactivity of PEEK. Int. J. Mol. Sci.; 2014; 15, pp. 5426-5445. [DOI: https://dx.doi.org/10.3390/ijms15045426] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24686515]

3. Panayotov, I.V.; Orti, V.; Cuisinier, F.; Yachouh, J. Polyetheretherketone (PEEK) for medical applications. J. Mater. Sci. Mater. Med.; 2016; 27, pp. 1-11. [DOI: https://dx.doi.org/10.1007/s10856-016-5731-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27259708]

4. Arif, M.F.; Kumar, S.; Varadarajan, K.M.; Cantwell, W.J. Performance of biocompatible PEEK processed by fused deposition additive manufacturing. Mater. Des.; 2018; 146, pp. 249-259. [DOI: https://dx.doi.org/10.1016/j.matdes.2018.03.015]

5. Kirin, B.S.; Kuznetsova, K.R.; Petrova, G.N.; Sorokin, A.E. Comparative analysis of the properties of polyetheretherketones of domestic and foreign production. Proc. VIAM; 2018; 5, pp. 34-43. [DOI: https://dx.doi.org/10.18577/2307-6046-2018-0-5-34-43]

6. Gurenkov, V.M.; Gorshkov, V.O.; Chebotarev, V.P.; Prudskova, T.N.; Andreeva, T.I. Polyetheretherketone. Properties, application, Production. Proceedings of the IV All-Russian Scientific and Technical Conference Polymer Composite Materials and Production Technologies of a New Generation; Moscow, Russia, 18 October 2019; pp. 51-66.

7. Lawrence, W.E.; Seferis, J.C.; Gillespie, J.W. Material response of a semicrystalline thermoplastic polymer and composite in relation to process cooling history. Polym. Compos.; 1992; 13, pp. 86-96. [DOI: https://dx.doi.org/10.1002/pc.750130204]

8. Burya, A.I.; Arlamova, N.T.; Chaika, L.V.; Hong, C. Application of kinetic equations to describe the thermal degradation of polyetheretherketone and carbon plastics based on it. Vestnik NovGU; 2015; 6, pp. 43-47.

9. Trufanova, N.M.; Subbotin, E.V.; Shcherbinin, A.G.; Kazakov, A.V.; Ershov, S.V. Mathematical modeling of heat and mass transfer processes during extrusion of nonlinear polymer compositions and experimental study of their properties. Bull. Perm Fed. Res. Cent.; 2017; 3, pp. 80-85.

10. PEEK Material Model. Available online: https://polymerfem.com/best-peek-material-model-2022/ (accessed on 1 November 2022).

11. Garcia-Gonzalez, D.; Rusinek, A.; Jankowiak, T.; Arias, A. Mechanical impact behavior of polyether–ether–ketone (PEEK). Compos. Struct.; 2015; 124, pp. 88-99. [DOI: https://dx.doi.org/10.1016/j.compstruct.2014.12.061]

12. Rae, P.J.; Brown, E.N.; Orler, E.B. The mechanical properties of poly(ether-ether-ketone) (PEEK) with emphasis on the large compressive strain response. Polymer; 2007; 48, pp. 598-615. [DOI: https://dx.doi.org/10.1016/j.polymer.2006.11.032]

13. Chen, F.; Ou, H.; Lu, B.; Long, H. A constitutive model of polyether-ether-ketone (PEEK). J. Mech. Behav. Biomed. Mater.; 2016; 53, pp. 427-433. [DOI: https://dx.doi.org/10.1016/j.jmbbm.2015.08.037] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26409233]

14. Qin, X.; Wu, X.; Li, H.; Li, S.; Zhang, S.; Jin, Y. Numerical and experimental investigation of orthogonal cutting of carbon fiber-reinforced polyetheretherketone (CF/PEEK). Int. J. Adv. Manuf. Technol.; 2022; 119, pp. 1003-1017. [DOI: https://dx.doi.org/10.1007/s00170-021-08317-3]

15. Chen, K.; Lu, Z.; Wei, P.; Liu, H.; Wei, D.; Xie, H. Study of Friction and Wear Characteristics of PEEK by Reciprocating Sliding Experiment and Temperature Dependences Simulation. Tribol. Lett.; 2022; 70, 99. [DOI: https://dx.doi.org/10.1007/s11249-022-01641-z]

16. Velisaris, C.N.; Seferis, J.C. Crystallization kinetics of polyetheretherketone (peek) matrices. Polym. Eng. Sci.; 1986; 26, pp. 1574-1581. [DOI: https://dx.doi.org/10.1002/pen.760262208]

17. Velisaris, C.N.; Seferis, J.C. Heat transfer effects on the processing-structure relationships of polyetheretherketone (PEEK) based composites. Polym. Eng. Sci.; 1988; 28, pp. 583-591. [DOI: https://dx.doi.org/10.1002/pen.760280907]

18. Barnes, J.A.; Simms, I.J.; Farrow, G.J.; Jackson, D.; Wostenholm, G.; Yates, B. Thermal expansion characteristics of PEEK composites. J. Mater. Sci.; 1991; 26, pp. 2259-2271. [DOI: https://dx.doi.org/10.1007/BF00549197]

19. Chapman, T.J.; Gillespie, J.W., Jr.; Pipes, R.B.; Manson, J.A.; Seferis, J.C. Prediction of Process-Induced Residual Stresses in Thermoplastic Composites. J. Compos. Mater.; 1990; 24, pp. 616-643. [DOI: https://dx.doi.org/10.1177/002199839002400603]

20. Turner, B.N.; Strong, R.; Gold, S.A. A review of melt extrusion additive manufacturing processes: I. Process design and modeling. Rapid Prototyp. J.; 2014; 20, pp. 192-204. [DOI: https://dx.doi.org/10.1108/RPJ-01-2013-0012]

21. Turner, B.N.; Gold, S.A. A review of melt extrusion additive manufacturing processes: II. Materials, dimensional accuracy, and surface roughness. Rapid Prototyp. J.; 2015; 21, pp. 250-261. [DOI: https://dx.doi.org/10.1108/RPJ-02-2013-0017]

22. Cao, Z.; Dong, M.; Liu, K.; Fu, H. Temperature Field in the Heat Transfer Process of PEEK Thermoplastic Composite Fiber Placement. Materials; 2020; 13, 4417. [DOI: https://dx.doi.org/10.3390/ma13194417]

23. Athale, M.; Park, T.; Hahnlen, R.; Pourboghrat, F. Experimental characterization and finite element simulation of FDM 3D printed polymer composite tooling for sheet metal stamping. Int. J. Adv. Manuf. Technol.; 2022; 121, pp. 6973-6989. [DOI: https://dx.doi.org/10.1007/s00170-022-09801-0]

24. Agelet de Saracibar, C.; Lundbäck, A.; Chiumenti, M.; Cervera, M. Shaped Metal Deposition Processes. Encyclopedia of Thermal Stresses; Springer: Dordrecht, The Netherlands, 2014; pp. 4346-4355. [DOI: https://dx.doi.org/10.1007/978-94-007-2739-7_808]

25. Lundbäck, A. Modelling of metal deposition. Finite Elem. Anal. Des.; 2011; 47, pp. 1169-1177. [DOI: https://dx.doi.org/10.1016/j.finel.2011.05.005]

26. Wang, P.; Zou, B.; Ding, S.; Li, L.; Huang, C. Effects of FDM-3D printing parameters on mechanical properties and microstructure of CF/PEEK and GF/PEEK. Chin. J. Aeronaut.; 2021; 34, pp. 236-246. [DOI: https://dx.doi.org/10.1016/j.cja.2020.05.040]

27. Zanjanijam, A.R.; Major, I.; Lyons, J.G.; Lafont, U.; Devine, D.M. Fused Filament Fabrication of PEEK: A Review of Process-Structure-Property Relationships. Polymers; 2020; 12, 1665. [DOI: https://dx.doi.org/10.3390/polym12081665]

28. Ramesh, M.; Rajeshkumar, L.; Balaji, D. Influence of Process Parameters on the Properties of Additively Manufactured Fiber-Reinforced Polymer Composite Materials: A Review. J. Mater. Eng. Perform.; 2021; 30, pp. 4792-4807. [DOI: https://dx.doi.org/10.1007/s11665-021-05832-y]

29. Belyaev, N.M.; Ryadno, A.A. Methods of Non-Stationary Heat Conduction: Proc. Allowance for Universities; 1st ed. M.: Higher school: Moscow, Russia, 1978; 328.

30. Pobedrya, B.E. Numerical Methods in the Theory of Elasticity and Plasticity: Proc. Allowance; 2nd ed. Publishing House of Moscow State University: Moscow, Russia, 1995; 366.

31. Standartinform. State All-Union Standard 4648-2014. PLASTICS. Three-Point Bending Flexural Test Method; Standartinform: Moscow, Russia, 2014.

32. Oskolkov, A.A.; Bezukladnikov, I.I.; Trushnikov, D.N. Rapid Temperature Control in Melt Extrusion Additive Manufacturing Using Induction Heated Lightweight Nozzle. Appl. Sci.; 2022; 12, 8064. [DOI: https://dx.doi.org/10.3390/app12168064]

33. Technical Data Sheet: CreatBot PEEK//CreatBot. 2022; 1p Available online: https://www.creatbot.com/downloads/filament/CreatBot%20PEEK-TDS%20EN.pdf (accessed on 1 November 2022).