# Determination of Crystallization Kinetics of Polyetheretherketone in Material Extrusion Process Using Numerical Simulation

Shahriar Bakrani Balani, Hossein Mokhtarian, Arthur Cantarel, France

Chabert, Eric Coatanéa, Valerie Nassiet

► To cite this version:

Shahriar Bakrani Balani, Hossein Mokhtarian, Arthur Cantarel, France Chabert, Eric Coatanéa, et al. Determination of Crystallization Kinetics of Polyetheretherketone in Material Extrusion Process Using Numerical Simulation. Polymers for Advanced Technologies, 2025, 36 (12), pp.e70458. 10.1002/pat.70458. hal-05489346

HAL Id: hal-05489346

https://ut3-toulouseinp.hal.science/hal-05489346v1

Submitted on 2 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

OPEN ACCESS

# Determination of Crystallization Kinetics of Polyetheretherketone in Material Extrusion Process Using Numerical Simulation

Shahriar Bakrani Balani<fnref n="1" /> ID | Hossein Mokhtarian<fnref n="1" /> ID | Arthur Cantarel² ID | France Chabert³ ID | Eric Coatanéa<fnref n="1" /> ID |
Valerie Nassiet³ ID

**Correspondence:** Shahriar Bakrani Balani (shahriar.bakranibalani@tuni.fi) | Hossein Mokhtarian (hossein.mokhtarian@tuni.fi)

**Received:** 30 September 2025 | **Revised:** 30 September 2025 | **Accepted:** 3 December 2025

**Keywords:** crystallization kinetics | fused filament fabrication | material extrusion | numerical simulation | PEEK

## ABSTRACT

Despite considerable progress over the past 3 years, some issues remain in obtaining printed polyetheretherketone (PEEK) parts printed by the material extrusion process (MEX). The fast crystallization of PEEK during printing results in distortion, poor dimensional accuracy, and low inter-filament adhesion. Optimizing the process requires a fine understanding of the material deposition and optimization of the features affecting the mechanical properties of the output. In this study, we modeled the polymer melt deposition on the substrate, heat transfer, and crystallization kinetics of PEEK. The isothermal and non-isothermal crystallization kinetics were determined for the processing range from the glass transition to the melting temperature. The non-isothermal crystallization of the polymer was applied to the process using the modified Nakamura non-isothermal crystallization equation. Moreover, the influence of the environmental temperature on the crystallization kinetics was determined for the deposition of the first layer on the substrate. This study determined the optimal process parameters to improve the mechanical properties of the manufactured parts, owing to a better interdiffusion of the polymer chains between the deposited layers and less porosity.

## 1 | Introduction

Printing polyetheretherketone (PEEK) is particularly challenging, due to its high processing temperature (613 K). Large internal stress, warpage, and delamination originate from a large temperature change during forming [1]. Arif et al. noted that PEEK specimens built in the vertical orientation exhibited delamination and poor mechanical performance, due to the high thermal gradient along the build direction [2].

PEEK and polyetherketoneketone (PEKK) are among the most promising thermoplastics used in the material extrusion process

(MEX) because their mechanical strength, chemical resistance, and thermal stability are suitable for applications in high-demanding applications such as the aerospace and medical fields. These remarkable properties mainly stem from the aromatic rings in the chemical structures and the crystallinity. Specifically, in the crystalline phase, the macromolecules are closely packed into spherulites, resulting in densification and improved thermal resistance of the entire material. PEEK exhibits a degree of crystallinity of approximately 40%, depending on the thermal history [3].

In addition to hobbyist 3D printing [4], the MEX process is a promising technology for industrial use. A solid filament enters

<footnote n="1">¹Faculty of Engineering and Natural Sciences, Tampere University, Tampere, Finland | ²Institut Clément (ICA), Université de Toulouse, CNRS UMR 5312, INSA, ISAE-Supaéro, INSA, IMT Mines Albi, UPS, Toulouse, France | ³Laboratoire Génie de Production, École Nationale D'ingénieurs de Tarbes (ENIT), University of Technology of Tarbes (UTTOP), Université de Toulouse, Tarbes, France</footnote>
a liquefier where it melts because of the heating electric resistance [5]. Then, the melted polymeric filament is deposited side-by-side and layer-by-layer to build the final part. Since 2015, PEEK has been considered a raw material suitable for the MEX process [6].

Most studies on PEEK printing by MEX process are driven by medical applications, which exploit the material biocompatibility and high mechanical strength. Additive manufacturing discloses new possibilities for manufacturing patient-specific implants, as manufactured parts can be readily available for use within a few hours after admission to the hospital. PEEK plates for general trauma and orthopedic surgeries have been successfully printed [7], as well as implants in cranio-maxillofacial surgery [8], dental and cranio-maxillofacial implants [9], lumbar spinal cages [10] and cranial implants [11]. Nevertheless, a recent review highlighted the lack of consensus on a suitable printing parameter combination for PEEK [12]. Optimized parameters for printing are therefore worth pursuing. Particularly, despite the extensive research on extrusion-based additive manufacturing, only a few studies dealt with the additive manufacturing of small PEEK parts, especially samples without defects such as warpage or delamination. Furthermore, PEEK properties are suitable for the harsh operating conditions required for spacecraft and aircraft, such as temperature changes and ultraviolet radiation. Composite-based PEEK is employed to print heat shields to protect spacecraft [13].

Many experimental studies on PEEK printing indicate that temperature is the dominant factor affecting print quality. Tseng et al. [14] revealed that the thermal conditions are the primary factor influencing the mechanical strength of PEEK specimens printed via MEX. Minimizing the thermal gradients across the beads can lead to PEEK parts with excellent macroscopic properties [2]. Thus, the nozzle temperature, substrate temperature, environment temperature, and cooling conditions during PEEK printing clearly influence crystallization and the final properties of the printed parts. Similar to the semicrystalline polymers, the properties of the final parts of PEEK strongly depend on the degree of crystallinity. The tensile strength of PEEK is 76 and 97 MPa at 16% and 40% degree of crystallinity, respectively [15].

Moreover, crystallization kinetics influence the auto-adhesion phenomenon, according to which each strand is tightly bonded with the adjacent strand through filament deposition. Auto-adhesion, also called coalescence, is based on the diffusion of macromolecules across the interface of two adjacent filaments. The fast crystallization hampers the macromolecule diffusion into the adjacent part, with a resulting poor bonding strength of PEEK [16]. Conversely, slow crystallization allows cross-interfacial crystallization, namely the link of spherulites located at the interface with macromolecules that were initially on each side of the interface. In turn, cross-interfacial crystallization promotes PEEK interfacial resistance, as reported by Cho and Kardos [17].

Yang et al. [18] investigated the correlation of the thermal processing conditions to crystallinity and mechanical properties in the MEX process. The environmental temperature changes in the range 300–475 K led to different degrees of crystallinity, even in different regions of the same PEEK part. The highest

percentage of crystalline phase was obtained at 473 K, along with a remarkable improvement in the elastic modulus and tensile strength. Moreover, when printing PEEK, a higher environment temperature assists in the slow crystallization and reduces internal stress. Wu proved that the warping deformation is minimal with a chamber temperature of 403 K and a nozzle temperature of 623 K [19]. To improve the printing conditions of PEEK, Hu added a heat collector and a high-power heater to the nozzle module and designed a nearly fully enclosed chamber [1]. The temperature profiles within 50 mm around the nozzle were approximately 323 K higher with the heat collector compared to printing in room temperature. As the sample printed with the traditional nozzle shows severe warpage (20.4%), the new printer design reduces the warpage rate to 5%. Higher crystallinity and no delamination of the PEEK parts were obtained, resulting in improved tensile and flexural properties. PEEK is characterized by higher crystallinity and better tensile properties when it is printed at higher nozzle temperatures, which allow for the appropriate polymeric chain crystallization [11, 18, 20]. For the samples processed at 693 K by Ding et al., the interlaminar bonding reached its maximum, which was attributed to the partial melting of the previous layer upon the molten filament deposition, thereby increasing the adhesion between layers [21]. However, some studies highlight the risk of degradation at high temperatures: Vaezi et al. observed a color change of PEEK when the nozzle temperature was in the range 703–723 K [22]. When PEEK is kept above 613 K under oxidative conditions, chain cleavage occurs adjacent to the carbonyl functional groups, and subsequent hydrogen abstraction from the aromatic ring leads to crosslinking between the adjacent aryl rings [23]. This recombination leads to a color change, a significant increase in shear viscosity, and a decrease in the material's crystallizability.

In summary, the nozzle temperature and environment temperature are decisive in controlling the crystallization and interfilament bonding [24]. During printing, the thermal gradient is more pronounced for PEEK than for other thermoplastics because of the high melting temperature. The thermal gradient causes weak interfacial adhesion, delamination, and warpage, leading to limitations in mechanical performance. Although the correlation between processing temperatures and mechanical properties is extensively investigated, insights into local phenomena are scanty. Particularly, detailed information on the crystallization occurring inside a single PEEK filament allows the fine-tuning of the crystallization and interfacial bonding within the whole part.

The thermal profile affects both the bead coalescence and the crystallization of the polymer upon cooling. Measuring heat transfer during the process is difficult because it requires non-destructive techniques. Seppala and Migler studied the heat transfer with acrylonitrile butadiene styrene (ABS) filaments using infrared thermography imaging [25]. The cooling rate depended on the printing speed and the number of layers; for the third layer, a printing speed of 10 mm s⁻¹ corresponded to a cooling rate of 60 K s⁻¹, whereas at a printing speed of 90 mm s⁻¹ the cooling rate was 180 K s⁻¹ for the same layer. For the eighth layer, the cooling rates at the same speeds were 50 and 10 K s⁻¹, respectively [25]. The combination of in situ infrared thermography and Raman spectroscopy was proposed by Northcutt in
2018 to measure the temperature and degree of crystallinity of extruded Polycaprolactam during deposition. Costanzo et al. [26] examined the thermal profile of the first layer using a micro-thermocouple placed on the printing plate in correspondence with the deposition area.

The challenge of determining the degree of crystallinity can be addressed through a modeling approach tailored for process optimization. The determination of the printed part crystallinity by numerical simulation is essential to optimize the reliability of the MEX process. The numerical simulation of the process requires fine modeling of the material deposition coupled with heat transfer. However, material deposition modeling is a primary challenge. While several studies have been conducted using computational fluid dynamics (CFD), where one-phase flow is coupled with the heat transfer equations to model the temperature variation in the liquefier [27, 28], a limited number of studies dealt with modeling the material deposition on the substrate in the MEX process. This is due to the difficulty in tracking the boundary of the polymer melt and air. Studies on the modeling of free-form extruders were carried out by Bonito et al. [29]. Viscoelastic flows were modeled with complex free surfaces in three space dimensions with a mathematical formulation similar to that of the volume of fluid (VOF) method [29]. The ability of the VOF method to track the polymer flow makes it suitable for modeling the polymer extrusion [30, 31]. Using VOF equations, Jun et al. accurately predicted the area of contact between adjoining strands for different speed ratios and layer thicknesses [32]. Dabiri et al. studied the deposition of hot viscous filaments with a model tracking the strand interfaces, but excluding the influence of the printing parameters on the flow behavior [33]. In 2015, Comminal modeled the behavior of the extrudate of a viscoelastic material with a numerical approach considering the viscoelasticity of the polymer and the streamlines owing to the elastic instabilities in the die [34]. The sharkskin phenomenon is related to the die elastic instability [27]. Comminal et al. modeled the extrusion of a strand of a semi-molten material on a moving plate [35]. The numerical simulation using CFD equations fully resolved the deposition flow and provided a cross-section of the printed strand. Furthermore, it was suitable to determine the influence of the gap distance and velocity ratio on the size and shape of the strand. This approach properly modeled the material deposition along a toolpath with a sharp corner and over-deposition of the material [36], enhancing the effect of toolpath smoothing on the shape of the bead-deposited material. According to the obtained results, a constant extrusion rate on the corners led to material over-deposition. Using an adaptable finite element analysis to determine the heat transfer in the MEX process, Amico et al. approached the material deposition as a uniform change in height across the part, instead of directly simulating the material flow [37]. The influence of the feed rate on the cooling rate of the deposited strand was included in their study. Furthermore, the thermal behavior of a RepRap 3D printer liquefier was studied by Jerez-Mesa et al., who modeled the fluid flow and temperature in the liquefier using continuity equations [38].

Level-set (LS) equations have gained attention for modeling multiphase systems [39]. The ability to track the fluid interface while considering the viscoelastic behavior of the fluid

is the main advantage of the LS equations. Therefore, the LS method is suitable for modeling free-form extrusion and studying coalescence, for example [24]. In our previous studies, we used the LS equations to model the free-form extrusion by considering the viscoelastic behavior of the polymer [39, 40]. Liravi et al. modeled the extrusion process of a fluid with low viscosity [41].

Despite the numerous experimental investigations determining the crystallinity of the printed parts, numerical simulation studies of crystallization in the MEX process are scarce. However, Li et al. predicted the crystallization and crystallization-induced warpage in an injection-molded semicrystalline polymer [42]. Anass et al. used the Schneider approach to model the crystallization coupled with a multiscale simulation by coupling the Navier-Stokes equations and LS. They described the nucleation and growth of crystalline entities under microinjection molding for a crystalline polymer [43]. Spina et al. proposed a multiscale and multi-physics simulation for semicrystalline thermoplastics [44], predicting the final microstructure of parts manufactured by injection molding under different processing conditions. McIroy and Graham proposed a modeling approach for flow-enhanced crystallization during the manufacturing of semicrystalline polymers with the MEX process [45].

To the best of our knowledge, the numerical modeling of in situ crystallization of PEEK in the MEX process is not reported in the literature. Due to the rapid crystallization of PEEK, determining its crystallinity with in situ measurement techniques is challenging. In this paper, a numerical model for depositing a high-viscosity polymer on a moving substrate is proposed. This model can be used to determine the influence of the printing parameters on the polymer flow, temperature field, and crystallinity. The outline of this paper is as follows (Figure 1): (i) prediction of the isothermal and non-isothermal crystallization of PEEK by applying the Hoffman-Lauritzen equation to the experimental results of crystallization determined by differential scanning calorimetry (DSC) and determination of the parameters of Avrami and Nakamura equations; (ii) numerical modeling of strand deposition on the moving substrate by applying TPF-LS equations; (iii) numerical determination of heat transfer and crystallinity of the deposited strand. The methodology developed in this study for PEEK is generalizable to other polymers and manufacturing processes.

## 2 | Materials and Methods

### 2.1 | Material Preparation

For this study, PEEK 450G in pellets provided by VICTREX was used. For all experiments, the polymer was dried for 24 h at 395 K in an oven to eliminate the moisture in the polymer structure [46]. After 24 h, the pellets were cooled to 300 K inside the furnace, and the pellets were kept in a desiccator to avoid water absorption.

The rheological, thermal, and thermomechanical analysis (TMA) samples were prepared via hot compression molding. The PEEK pellets were placed in a mold and heated above the melting temperature (633 K) at 70 kN pressure. The temperature
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Flowchart of the procedure of this study to determine the crystallization kinetics in the MEX process.</figcaption></figure>

and pressure were maintained for 300 s. Then, the mold was cooled under the same pressure at a cooling rate of 10 K min⁻¹.

## 2.2 | Characterization Techniques

The DSC tests were performed with a Q200-TA instrument under nitrogen flow. All experiments were carried out at different isothermal temperatures between 583 and 595 K, with an approximate sample weight of 10 mg. The samples were heated above the melting point to remove their prior thermomechanical history. This step was also essential for simulating the conditions inside the liquefier during the MEX process. The polymer was fully melted in the liquefier before being deposited layer-by-layer during the fabrication process. The glass transition (T_g), melting temperature (T_m), and enthalpy of crystallization (ΔH_c) were measured for each isothermal temperature.

The density and its variation with temperature were measured by TMA using a Perkin Elmer TMA7 in the temperature range 300–520 K with a temperature ramp of 3 K min⁻¹ under Helium flow. Samples with dimensions of approximately 5 × 5 × 2 mm were used. The tests were carried out by applying a constant force of 10 mN using a 3.7 mm diameter tip quartz probe. The thermal conductivity and thermal diffusivity of PEEK were determined using a TPS 2500 hot-disk thermal analyzer. Three measurements of each sample were performed to ensure reproducibility. The rheological characterization was carried out using an ARES LN2 rheometer equipped with an electrical oven under air flow conditions using a 25 mm diameter parallel-plate configuration. First, time sweep tests were performed to verify PEEK's thermal stability within the duration of the experiments (approximately 20 min). Strain sweeps were then applied at 1 rad s⁻¹ to delimit the viscoelastic linear domain, and frequency sweeps were performed at 1% within the range 100–0.05 rad s⁻¹.

## 3 | Experimental Analysis

The proposed simulation in this study requires the inclusion of certain properties of the PEEK material and the air environment, namely, thermal properties, density, and rheological properties. PEEK properties were experimentally measured, while the corresponding data for air were obtained from the literature. A thorough characterization of the isothermal and non-isothermal crystallization of PEEK is presented in this section.

### 3.1 | Determination of the Thermal Properties

Using a hot-disk analysis, the temperature variation of the samples was measured during the energy transfer from the transient plane source over a specific time interval. The thermal conductivity (λ) of PEEK was determined from 294 to 493 K. The thermal diffusivity was further calculated considering the density variation and thermal capacity variation with temperature. Figure 2 represents the PEEK thermal conductivity and thermal diffusivity in the range of 294 to 493 K. Thermal conductivity and thermal diffusivity have been extrapolated for temperatures above 500 K and melting temperature.

As shown in Figure 2a, PEEK thermal conductivity at 300 K was 0.286 W m⁻¹ K⁻¹ and it increased with temperature to reach a plateau at 0.330 W m⁻¹ K⁻¹ near 425 K, corresponding to the glass transition. At this temperature, the molecular motion is amplified, causing heat loss in all directions within the sample. The measurement was not possible above 500 K with the used apparatus.

Thermal diffusivity strongly depends on the thermal energy distribution through the material, consequently on the material homogeneity. The thermal diffusivity was approximately 0.22 mm² s⁻¹ (Figure 2b). Unlike thermal conductivity, the
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Correlation of the thermal properties with temperature. (a) Thermal conductivity and (b) thermal diffusivity.</figcaption></figure>

<table><caption>**TABLE 1** | Properties of air at 294 K.</caption><thead><tr><th>Parameters</th><th>Unit</th><th>Air</th></tr></thead><tbody><tr><td>Thermal conductivity</td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>0.024</td></tr><tr><td>Density</td><td>kg m<sup>-3</sup></td><td>1.225</td></tr><tr><td>Heat capacity</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>1005</td></tr><tr><td>Thermal diffusivity</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>1.95 &times; 10<sup>-5</sup></td></tr><tr><td>Viscosity</td><td>Pa s</td><td>1.810<sup>-6</sup></td></tr></tbody></table>

thermal diffusivity decreased at higher temperatures. The thermal diffusivity was calculated from the density and heat capacity variations. Therefore, the thermal diffusivity changes were caused by the great variation in density, while the polymer went through its glass transition. The activation energy and density of PEEK are 82.8 kJ mol<sup>-1</sup> and 1320 kg m<sup>-3</sup>, respectively [46, 47].

The thermal properties, density, and viscosity of air at room 300 K, obtained from the literature, are listed in Table 1 [48, 49].

## 3.2 | Shear-Dependent Viscosity

Rheological properties of PEEK have been studied in our previous studies [50]. The shear-dependent Carreau-Yasuda model presented in Equation (1) has been used for the viscosity in the numerical model.

$$ \eta = \eta_{\text{inf}} + (\eta_0 - \eta_{\text{inf}}) [1 + (\lambda\dot{\gamma})^a]^{\frac{n-1}{a}} \quad (1) $$

where $\eta_0$ is the viscosity of the fluid at zero shear rate, $\eta_{\text{inf}}$ is the viscosity of the fluid at infinite shear rate, $\lambda$ is the relaxation time, $n$ is the power index, $a$ is a dimensionless parameter, $n$ is the pseudoplasticity index, and $\dot{\gamma}$ is the shear rate. The terms of the Carreau-Yasuda model of PEEK used in the numerical simulations are presented in Table 2.

<table><caption>**TABLE 2** | Terms of Carreau model for PEEK at 656 K.</caption><thead><tr><th>Polymer</th><th>PEEK</th></tr></thead><tbody><tr><td>η<sub>0</sub></td><td>7071 ± 153</td></tr><tr><td>η<sub>inf</sub></td><td>0</td></tr><tr><td>λ</td><td>1.45 ± 0.6</td></tr><tr><td>a</td><td>0.78 ± 0.12</td></tr><tr><td>n</td><td>0.59 ± 0.06</td></tr></tbody></table>

## 3.3 | Isothermal Crystallization

Isothermal crystallization was studied prior to PEEK non-isothermal crystallization in the MEX process. The degree of crystallinity of a semicrystalline polymer is defined as the ratio of the endothermic crystallization enthalpy to the crystallization enthalpy of the same polymer at its fully crystalline state (Equation 2). The enthalpy of the fully crystalline PEEK is $\Delta H_{100\%} = 130$ J g<sup>-1</sup> [51]. The enthalpy of crystallization is the area under the crystallization peak determined by the DSC tests.

$$ \alpha(t) = \frac{\int_{0}^{t} \delta h}{\Delta H_{100}} \quad (2) $$

where $\alpha(t)$ is the degree of crystallinity at different times during isothermal crystallization, $\delta h$ is the crystallization enthalpy, and $\Delta H_{100}$ is the enthalpy for the full crystallization of the polymer.

The kinetics of the crystallization of polymers under isothermal conditions is well described by the Avrami equation [52, 53] (Equation 3), which determines the correlation of the degree of crystallinity at an isothermal temperature and time. At each isothermal temperature, the parameters of the Avrami equation, namely the coefficient $k_{av}$ and the Avrami exponent $n$, were determined separately.
$$ \alpha(t) = 1 - \exp(-k_{av}t^n) \quad (3) $$

From Equation (4), which is the logarithmic form of Equation (3), $\ln(-\ln[1 - \alpha(t)])$ was plotted versus $\ln(t)$ for each isothermal temperature (Figure 8) with $n$ as the slope of the curve and $k_{av}$ the intercept with the ordinate axis.

$$ \ln(-\ln[1 - \alpha(t)]) = \ln k_{av} + n \ln t \quad (4) $$

The terms of the Avrami equation can be obtained by analyzing the DSC test results related to the temperature range from $T_g$ to $T_m$. However, for PEEK, the crystallization kinetics can be extremely fast, and the isothermal crystallization for certain temperatures exceeds the measurement limits of the DSC apparatus used in this study. Therefore, even at very fast cooling rates, from the equilibrium melting state to the targeted isothermal temperature, the polymer crystallizes during the cooling ramp. Consequently, the determination of the Avrami terms by DSC for these polymers is only possible when the crystallization kinetics are relatively slow. Experimental results showed very slow crystallization kinetics near $T_g$ and $T_m$. Figure 3 shows the peak related to PEEK non-isothermal crystallization (cooling ramp) obtained by DSC at a cooling rate of 30 K min$^{-1}$, at which the polymer crystallization starts at 583 K and finishes at 524 K, and is completely absent above 583 K. Therefore, in the case of PEEK 450G, we could isothermally study the crystallization for temperatures above 583 K by DSC. However, to determine the non-isothermal crystallization in the MEX process, the Avrami terms for all the temperature ranges from $T_g$ until $T_m$ should be obtained.

The temperature dependence of the Avrami kinetics of the crystallization constant ($k_{av}(T)$) was modeled using the Hoffman-Lauritzen equation (Equation 5), according to which the crystallization is symmetrical on both sides of the half-time crystallization ($t_{1/2}$).

$$ k_{av}(T) = \left(\frac{1}{t_{1/2}}\right) = K_0 \exp\left(\frac{-U}{R(T - T_{\infty})}\right) \exp\left(\frac{-K_G(T + T_0)}{2T^2 \Delta T}\right) \quad (5) $$

<figure><img src="image_5.png" /><figcaption>FIGURE 3 | Crystallization of PEEK 450G with 30 K min$^{-1}$ cooling rates. Crystallization starts at $T = 583$ K.</figcaption></figure>

The first exponential term refers to the driving force of crystallization, whereas the second exponential term refers to the molecular chain mobility. $U$ is the diffusional activation energy for the transport of crystallizable segments at the liquid-solid interface [54], $R$ is the universal gas constant, and $K_G$ is the nucleation characteristic that influences the position of the axial peak of the $k_{av} \cdot T_{\infty}$ is the temperature at which the motion associated with viscous flow ceases (typically 30 K below $T_g$) [55]. $T_0$ is the equilibrium melting temperature, which is obtained by the linear Hoffman-Week method for PEEK [56, 57]. $K_G$ and $K_0$ were determined by interpolating the Hoffman-Lauritzen equation from experimental data. The Hoffman-Lauritzen model was used to determine the Avrami terms below 583 K. Specifically, the Hoffman-Lauritzen equation extrapolation was used to predict PEEK crystallization in the entire crystallization range from $T_m$ until $T_g$ (Table 3).

## 3.4 | Non-Isothermal Crystallization

According to Nakamura et al. [58], the ratio of the transformed phase from amorphous to solid (semicrystalline) state for a semicrystalline polymer under non-isothermal conditions is obtained using Equation (6):

$$ \alpha(t) = 1 - \exp\left(-\left(\int_{0}^{t} K(T) d\tau\right)^{n}\right) \quad (6) $$

where $\alpha(t)$ is the fraction of the polymer converted into the semicrystalline state, $t$ is the time, $n$ is the Avrami exponent analogous to the one of the isothermal studies, and $K(T)$ is the Nakamura kinetics of the crystallization coefficient. Reformulation of Equation (6) results in Equation (7).

$$ \frac{d\alpha}{dt} = nK(T)(1 - \alpha)\left[\ln\left(\frac{1}{1 - \alpha}\right)\right]^{\frac{n-1}{n}} \quad (7) $$

Furthermore, $K(T)$ could be obtained from $k_{av}(T)$ or $t(T)^{\frac{1}{n}}$ according to Equation (8).

$$ K(T) = k_{av}(T)^{\frac{1}{n}} = \ln(2)^{\frac{1}{n}} * \frac{1}{t(T)^{\frac{1}{n}}} \quad (8) $$

Crystallization is an exothermic transition that leads to a temperature increase in the system. We integrated the released

<table><caption>TABLE 3 | Value of the constants in the Hoffman-Lauritzen equation for PEEK 450G.</caption><thead><tr><th>Constants</th><th>Unit</th><th>PEEK 450G</th></tr></thead><tbody><tr><td>Diffusional activation energy (U)</td><td>J mol<sup>-1</sup></td><td>6300</td></tr><tr><td>Universal gas constant (R)</td><td>J mol<sup>-1</sup> K<sup>-1</sup></td><td>8.314</td></tr><tr><td>(T<sub>∞</sub>)</td><td>K</td><td>400</td></tr><tr><td>Equilibrium melting temperature (T<sub>0</sub>)</td><td>K</td><td>663</td></tr></tbody></table>
heat due to crystallization by using PEEK specific heat capacity ($C_p$) and adding a $\Delta T_c$ in the heat equation, defined as follows:

$$ \Delta T_c = \frac{\Delta H}{C_p} * \frac{d\alpha}{dt} \qquad (9) $$

where $\Delta H$ is the crystallization enthalpy. In Equation (9), $C_p$ is considered constant and independent from temperature.

To consider the fluid flow in the system, a convection term was added to the Nakamura crystallization equation (Equation 7). Afterwards, we considered the fluid and transferred the crystalline phase during the fluid flow. Finally, the non-isothermal crystallization for a polymer fluid flow was determined using Equation (10):

$$ \frac{\partial \alpha}{\partial t} + (u \cdot \nabla) \alpha = nK(T)(1 - \alpha) \left[ \ln \left( \frac{1}{1 - \alpha} \right) \right]^{\frac{n-1}{n}} \qquad (10) $$

where $\alpha$ is the degree of relative crystallinity, and $u$ is the velocity field of the fluid.

## 3.5 | PEEK Crystallization Kinetics

The first step to determine the isothermal crystallization was selecting the measurable temperature range by DSC. Particularly, we determined the crystallization starting temperature because PEEK crystallization during the isothermal study must be complete under isothermal conditions. A non-isothermal test using a DSC apparatus was carried out. For the highest cooling rates, the crystallization temperature shifted toward lower temperatures (up to $T_g$) and at the lowest cooling rates, up to $T_m$. Consequently, a cooling rate of 2 K min⁻¹ was selected to determine the starting temperature for crystallization. During the cooling step, the crystallization started at approximately 583 ± 3 K (Figure 4). Accordingly, the isothermal test was performed at temperatures higher than 583 K. Nevertheless, at temperatures higher than 596 K, PEEK crystallization was slower. Therefore, the isothermal tests were performed between 583 and 596 K.

The isothermal crystallization of PEEK 450G was performed in the range of 595 to 581 K. The isothermal tests at $T = 588$ K shown in Figure 5, included four steps: the sample was initially heated at 648 K at a heating rate of 50 K min⁻¹ and kept at this temperature for 240 s to eliminate the polymer's prior crystalline history. Afterwards, the polymer was cooled down at a rate of 30 K min⁻¹ until the isothermal target temperature was reached. The polymer was kept at its isothermal temperature until crystallization was completed.

PEEK isothermal crystallization was performed at nine different isothermal temperatures (Figure 6). The results for 581, 583, and 585 K show that a small fraction of the amorphous phase undergoes non-isothermal crystallization. Consequently, these three isothermal temperatures were not considered in this study.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 | DSC curve of the PEEK 450G. The initial heating rate for eliminating the initial crystallinity was 10 K min⁻¹; samples were cooled and heated at 2 K min⁻¹ during the second heating cycle.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 5 | Isothermal DSC study of PEEK 450G at $T = 588$ K.</figcaption></figure>

Table 4 lists the enthalpy, degree of crystallinity, and half-time crystallization obtained from the isothermal studies. The variation of the crystallization enthalpy was between 25.3 J g⁻¹ equivalent to 19.5% and 31.1 J g⁻¹ equivalent to 24% of the crystallinity. The half-time crystallization increased as the temperature decreased, since polymer crystallization was slower at higher temperatures.

Below 585 K, the crystallization points of the Hoffman-Lauritzen equation were fitted on the half-time determined experimentally using Origin software, resulting in a good match (Figure 7). Table 5 reports the values and standard deviation of the parameters of the Hoffman-Lauritzen equation for PEEK obtained by fitting the curve to the experimental points.

After determining the Avrami kinetics by the half-time crystallization (Equation 8), Avrami kinetics of the crystallization
coefficient and Avrami exponent were determined for PEEK by applying the Avrami approach. For each isothermal temperature in the range of 595 to 588 K, the curve of $ln(-ln[1 - α(t)])$ versus $ln(t)$ was plotted in Figure 8. The fitting curve is a straight line with slope *n* and intercept $k_{av}$. The Avrami kinetics and exponents for the isothermal temperatures are listed in Table 6.

The Avrami parameters provide information on the crystalline structure of the polymer, especially the crystal formation during crystallization [54]. As shown in Table 6, the Avrami exponent (*n*) varied from 3.18 to 3.74, with a value approaching 3 at lower temperatures and 4 at a temperature higher than 595 K. When the value of *n* is approximately 3, the shape of the crystallization formation is 2D spherulites [51], whereas *n* = 4 indicates 3D spherulitic formation. At temperatures close to $T_m$ or $T_g$, the decrease in crystallization kinetics becomes more evident. As shown in Table 6, at 588 K, the kinetics of crystallization were faster than at 595 K (closer to $T_m$).

Figure 9 shows the good agreement of the Avrami kinetics coefficient determined by the Avrami approach and the predictive model of Equation (8). At a temperature close to $T_m$, the Avrami exponent approaches 4, whereas at lower temperatures it tends toward 3. Therefore, in the predictive model, we used the Avrami exponent from the experimental study for the range 595–588 K, and

<figure><img src="image_8.png" /><figcaption>FIGURE 6 | Isothermal study of PEEK 450G at different temperatures from 595 to 581 K.</figcaption></figure>

we assumed the Avrami exponents 4 and 3 for temperatures higher than 595 K and lower than 588 K, respectively. This assumption is consistent with our observations on the presence of two crystalline structures during the melting of PEEK. The DSC results showed that the transition temperature from the first to the second crystalline structure occurred at approximately 588 K (Figure 3). Figure 9 shows that the maximum value of the kinetics of crystallization for PEEK 450G is around 516 K. However, close to $T_m$ and $T_g$, the kinetics of crystallization highly decreases.

The Avrami kinetics of crystallization $k_{av}(T)$ and Avrami exponent *n* were inserted in Equation (8) to determine the Nakamura kinetics of crystallization $K(T)$ for the grade of PEEK we chose in this study. In Figure 10, we compare the Nakamura crystallization kinetics coefficients obtained for PEEK 450G to those found in the literature for another grade with lower viscosity, PEEK 150G, determined experimentally [59].

<figure><img src="image_9.png" /><figcaption>FIGURE 7 | Half-time crystallization and Hoffman-Lauritzen fitting of the experimental points.</figcaption></figure>

<table><caption>TABLE 5 | Parameters of the Hoffman-Lauritzen equation determined for PEEK.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Standard error</td></tr></thead><tbody><tr><td>K<sub>0</sub></td><td>7.68 × 10<sup>8</sup></td><td>1.09 × 10<sup>7</sup></td></tr><tr><td>K<sub>G</sub></td><td>-8.87 × 10<sup>6</sup></td><td>5.9 × 10<sup>5</sup></td></tr></tbody></table>

<table><caption>TABLE 4 | Influence of different isothermal studies on the enthalpy and half-time crystallization of PEEK 450G.</caption><thead><tr><td>Isothermal temperature (K)</td><td>Sample mass (mg)</td><td>Enthalpy ±0.2Jg<sup>-1</sup></td><td>Degree of crystallinity ±0.2%</td><td>Half-time crystallization (s)</td></tr></thead><tbody><tr><td>588</td><td>10.4</td><td>30.1</td><td>23.1</td><td>98</td></tr><tr><td>589</td><td>8.9</td><td>30.9</td><td>23.8</td><td>199</td></tr><tr><td>590</td><td>10.9</td><td>28.8</td><td>22.2</td><td>230</td></tr><tr><td>591</td><td>12.6</td><td>31.1</td><td>24.0</td><td>290</td></tr><tr><td>593</td><td>9.2</td><td>30.2</td><td>23.0</td><td>404</td></tr><tr><td>595</td><td>9.6</td><td>25.3</td><td>19.5</td><td>586</td></tr></tbody></table>
# 4 | Modeling and Numerical Simulation of Polymer Deposition in the MEX Process

## 4.1 | Flow Model

A Eulerian description of the flow field was considered. The 2D numerical simulation of the deposition process was carried out using a two-phase flow (TPF) simulation of COMSOL Multiphysics software. Mass conservation (Equation 11) and conservation of momentum with Navier-Stokes equations (Equation 12) has been used to simulate the flow, advection, and follow the front of the polymer [60, 61].

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho u) = 0 \qquad (11) $$

$$ \rho \left( \frac{\partial u}{\partial t} + u \cdot \nabla u \right) = \nabla \cdot \left[ -pI + \eta(T) \left( \nabla u + (\nabla u)^T - \frac{2}{3} (\nabla \cdot u) I \right) \right] + \rho g + F_{\text{st}} + F_e \qquad (12) $$

In these equations, $\rho$ is the density, $u$ is the flow velocity, $p$ is the pressure applied to the fluid, $\eta$ is the dynamic viscosity of the fluid, $g$ is the gravity field, $F_{\text{st}}$ represents the force resulting from surface tension, and $F_e$ indicates all the other external forces. The Navier-Stokes equations correspond to the contributions of different forces applied to the fluid.

For flows at very low Reynolds numbers (for polymer flow in the MEX nozzle, the Reynolds number is less than 0.1), the inertial term in the Navier-Stokes equation can be neglected. The Stokes equation used in this study is Equation (13).

<figure><img src="image_10.png" /><figcaption>**FIGURE 8** | Determination of the Avrami kinetics $k_{av}$ and exponents $n$ by the $\ln$ method.</figcaption></figure>

<table><caption>**TABLE 6** | Avrami isothermal kinetics and exponent for the different isothermal temperatures.</caption><tr><td>Temperature (K)</td><td>595</td><td>593</td><td>591</td><td>590</td><td>589</td><td>588</td></tr><tr><td>Avrami exponent $n$</td><td>3.55</td><td>3.63</td><td>3.74</td><td>3.49</td><td>3.35</td><td>3.18</td></tr><tr><td>Avrami kinetics $k_{av}$ (s<sup>-n</sup>)</td><td>1.46 × 10<sup>-11</sup></td><td>2.97 × 10<sup>-11</sup></td><td>6.75 × 10<sup>-11</sup></td><td>6.20 × 10<sup>-10</sup></td><td>2.06 × 10<sup>-9</sup></td><td>1.06 × 10<sup>-8</sup></td></tr></table>

$$ \rho \frac{\partial u}{\partial t} = \nabla \cdot [-pI + \eta(T)(\nabla u + (\nabla u)^T)] + \rho g + F_{\text{st}} + F_e \quad (13) $$

## 4.2 | Front Tracking

In the TPF simulation, the interface of the two fluids is tracked by the LS method [62], with the transport equation reported in Equation (14):

$$ \frac{\partial \varphi}{\partial t} + \nabla \cdot (u \varphi) = \gamma \nabla \cdot \left( \varepsilon_{\text{Is}} \nabla \varphi - \varphi(1-\varphi) \frac{\nabla \varphi}{|\nabla \varphi|} \right) \quad (14) $$

<figure><img src="image_11.png" /><figcaption>**FIGURE 9** | Determination of the Avrami kinetics along the processing range $T_g$ (433 K) to $T_m$ (623 K).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 10** | Nakamura kinetics of crystallization for PEEK 450G compared with PEEK 150G [59].</figcaption></figure>
where $u$ is the flow velocity; $\varphi$ is the volume fraction of the polymer, $\gamma$ is a re-initialization parameter, and $\varepsilon_{Is}$ is the parameter controlling the interfacial thickness. The re-initialization parameter $\gamma$ is considered as the maximum or close to the maximum velocity of the fluid in the TPF system to ensure the consistency of the results with the entire simulation. Lowering the $\varepsilon_{Is}$ value caused a reduced thickness of the interface to improve the clarity between the phases. The volume fraction was the output of the level-set method used to determine the position of the interface between the two phases and the properties of the meshes. The value of the color function $\varphi$ equal to 1 corresponds to a polymer in the system or in this mesh, whereas the value of the volume fraction equal to 0 indicates the presence of air in the system or in the mesh.

In the first step of the TPF simulation, the density and viscosity of each mesh of the system were determined using the mixing law in Equation (15), according to the volume fraction $\varphi$ obtained by the level-set method:

$$ \rho = \varphi\rho_2 + (1 - \varphi)\rho_1 \\ \eta = \varphi\eta_2 + (1 - \varphi)\eta_1 \qquad (15) $$

where $\rho_i$ and $\eta_i$ are the density and the viscosity of phase $i$, respectively.

## 4.3 | Heat Transfer

The partial differential heat equation was added to the model to reproduce the heat transfer in the flow simulation [63]. To determine the temperature variation in our model, the classical advection-diffusion in Equation (16) considered the fluid flow in the system and a source term $Q$ to include the contribution of crystallization or melting enthalpy:

$$ \frac{\partial T_s}{\partial t} + (u \cdot \nabla) T_s = \frac{\lambda}{\rho C_p} \nabla^2 T_s + Q \qquad (16) $$

where $T$ is the temperature, $u$ is the fluid velocity that links the heat transfer equation to the TPF simulation, $\lambda$ is the thermal conductivity, $\rho$ is the density, $C_p$ is the heat capacity of the fluid, and $Q$ is the heat flux connected to crystallization or melting.

## 4.4 | Boundary Conditions

The boundary conditions used for the numerical simulation of the material deposition on the substrate are presented in Figure 11.

To assign the proper thermal properties to each phase (polymer and air), three different methods were examined in the simulation to accurately determine the heat transfer in the MEX process:

* Mixing law (Equation 15);

* Step function (Equation 17);

* A combination of Equations (17) and (15). In this method, the thermal diffusivity was inserted to combine the two methods. The density and specific heat were introduced as the mixing law, and thermal conductivity was introduced in the model as a step function.

$$ \begin{array}{l} \text{If } \varphi \le 0.5 \text{ then properties of air else} \\ \text{If } 0.5 \le \varphi \le 1 \text{ properties of polymer} \end{array} \qquad (17) $$

The combination of mixing law and step function led to a correct phase assignment of the material properties [64].

# 5 | Results and Discussion

## 5.1 | Numerical Study of the Polymer Deposition in the MEX Process

This section reports the influence of the temperature of the outlet, substrate, and air on the heat transfer and degree of crystallinity of the first deposited strand of PEEK during the deposition phase of the MEX process. The determination of the optimal crystallization temperature is necessary to print PEEK parts with higher mechanical strength.

The determination of the temperature and degree of crystallization required the fine modeling of the material deposition. Consequently, prior to studying heat transfer, the material deposition and polymer flow were modeled. As shown in Figure 12, material deposition of the first strand started at $t=0$s and was

<figure><img src="image_13.png" /><figcaption>FIGURE 11 | Illustration of the boundary conditions and geometry used for numerical simulation.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**FIGURE 12** | Numerical simulation of the first layer deposition on the moving substrate at different moments.</figcaption></figure>

completed at 4.7s. The total length of the strands was 105 mm. The mean relative velocity of the substrate was $23 \text{ mm s}^{-1}$ for 4.7s.

To study the temperature influence on crystallinity, five different case studies were selected. Since the primary objective of this investigation was to determine the influence of temperature on the crystallization process, case studies 1, 2, and 3 were defined. The most significant and most common temperatures for PEEK, namely 303 and 433 K, were selected. Particularly, printing at room temperature is $T=303 \text{ K}$ and $T=433 \text{ K}$ is the glass transition temperature ($T_g$) of PEEK. Printing at room temperature is the case of material deposition in a chamber without any heating during printing, which is the typical scenario for home printing, simpler than controlling the build temperature. Contrarily, $T=433 \text{ K}$ is PEEK $T_g$, selected because the polymer is in the viscoelastic domain and molecular chains have mobility above this temperature.

Essential objective is to determine the optimal build temperature at which crystallization has not yet started at the end of the material deposition of the strand. To achieve the optimal mechanical strength, the following strand must be deposited when it is still in its amorphous state. However, the time between the end of the deposition of the first strand and the deposition of the second strand depends on the printing parameters, the deposition strategy, and the size of the part. Consequently, the build temperatures for the two case studies (case 4 and case 5) were determined in two stages (Table 7), crystallization after 4.7 and 12s, providing the appropriate temperature range for the optimal temperature.

Figure 12 shows the evolution of $\phi$, the volume fraction of the polymer in the system during the first layer deposition on the substrate. Determination of $\phi$ is prerequisite for modeling heat transfer and crystallization.

Figure 13 shows the temperature field during the material deposition for Case 3 at different times. During the material deposition, most of the extrudate's center, partially in the melted state, stayed above 565 K. The results obtained also show that the temperature at the front of the strand sharply decreases after exiting from the nozzle. Furthermore, the

<table><caption>**TABLE 7** | Case studies for the numerical simulation of temperature distribution and crystallization rate.</caption><thead><tr><td>Case studies</td><td><i>T</i><sub>m</sub> (K)</td><td><i>T</i><sub>sub</sub> (K)</td><td><i>T</i><sub>env</sub> (K)</td></tr></thead><tbody><tr><td>Case 1</td><td>661</td><td>303</td><td>303</td></tr><tr><td>Case 2</td><td>661</td><td>433</td><td>303</td></tr><tr><td>Case 3</td><td>661</td><td>433</td><td>433</td></tr><tr><td>Case 4</td><td>661</td><td>558</td><td>558</td></tr><tr><td>Case 5</td><td>661</td><td>573</td><td>573</td></tr></tbody></table>

exothermic crystallization reaction slowed the cooling of the strand.

Figure 14 shows the temperature variation at the centerline of the deposited strand at different times for Case 3. As time increased, the temperature decreased. However, the cooling rate at the center of the extrudate was relatively slow and promoted the crystallization of the polymer. As shown in Figure 14, crystallization influenced the temperature curves.

Figure 15 presents the evolution of the degree of crystallinity of the extrudate for Case 3. The amorphous phase and air are shown in dark blue. From light blue to orange (max 23%), different percentages of crystallinity were observed. At the extrudate front, where the cooling rate is high, crystallinity reaches its maximum (approximately 23%). However, at the center of the deposited bead, where the temperature is still high, the polymer is still in an amorphous state.

The degree of crystallinity field for the three studied cases at $t=12 \text{ s}$ is compared in Figure 16. The results show the influence of the environment and substrate temperatures on the degree of crystallinity of the deposited strand for the three studied cases. At $t=12 \text{ s}$, the crystallization of the strand was divided into three different regions.

Region 1 was the front part of the extrudate that rapidly cooled after exiting from the nozzle. The crystallization reached its maximum value early; therefore, for this part, the coalescence and interdiffusion of the deposited beads are disfavored.
<figure><img src="image_15.png" /><figcaption>FIGURE 13 | Temperature field for Case 3, namely $T_{sub} = T_{env} = 433$ K.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 14 | Temperature field at different moments on the extrudate centerline for Case 3.</figcaption></figure>

Region 2 underwent crystallization. The analysis of the kinetics of crystallization for the three studied cases showed that the crystallization rate for Case 1 is much faster than that for the other cases. The maximum values of the crystallization rate for Cases 1, 2, and 3 were $0.107\% \text{ s}^{-1}$, $0.215\% \text{ s}^{-1}$, and $0.325\% \text{ s}^{-1}$, respectively. At the bead surface, where the polymer in contact with air and substrate had the highest cooling rate, the crystallization was relatively fast and reached its maximum value before the crystallization of the center of the strand. In Region 3, the center of the extrudate was still at the melting temperature, and no crystallization was observed at 12s. The part of the bead just after exiting from the nozzle did not have enough time for cooling and crystallization. However, at the surface where the cooling rate was the highest, polymer crystallization was observed. For Case 1, the crystallization rate was between $10^{-3}\% \text{ s}^{-1}$ and $10^{-2}\% \text{ s}^{-1}$. However, for Case 2 and Case 3, the crystallization rate was negligible (less than $10^{-5}\% \text{ s}^{-1}$).

Figure 17 shows the temperature distribution and degree of crystallinity at the center of the deposited strand for the three cases at $t=12$ s. The influence of the crystallization on the temperature field was observed.

## 5.2 | Determination of the Optimal Environment Temperature to Promote the Interdiffusion of Strands and Layers

To improve the mechanical properties of the printed parts, the interdiffusion between two consecutive deposited layers was promoted and porosity prevented. The polymer bead remained in its amorphous state and above $T_g$. Therefore, the cooling rate of the deposited bead was controlled to preserve the amorphous phase and allow the interdiffusion before crystallization. This can be done using several methodologies: (i) Increasing the temperature of the chamber (or substrate) to reduce the cooling rate; (ii) Increasing the melting point in the liquefier (yet below the temperature at which the polymer degrades); (iii) Increasing the travel speed of the nozzle to quickly deposit the bead one upon the other before the cooling of the deposited bead; and (iv) Modifying the deposition strategy to drop the strands next to each other as fast as possible.

The travel speed and deposition strategy are independent of the material used and can be modified during printing. Conversely, the temperature of chamber (or substrate) and melting temperature are correlated to the material properties. Therefore, increasing melting temperature can be achieved by increasing the power of the liquefier heaters. Consequently, PEEK printers should operate at approximately 673 K.

Two critical durations were determined: the first bead deposition time of 4.7s and the second bead deposition time of 12s (go/back/go).

The numerical simulation showed that at a chamber temperature of 558 K, the crystallization of the surface of the strand
<figure><img src="image_17.png" /><figcaption>**FIGURE 15** | Distribution of the degree of crystallinity for Case 3—T<sub>sub</sub> = T<sub>env</sub> = 433 K.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 16** | Non-isothermal crystallization of PEEK at *t* = 12 s for the different case studies.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**FIGURE 17** | Comparison of temperature and degree of crystallinity in the center of the deposited extrudate for Cases 1, 2, and 3 after 12 s.</figcaption></figure>

started at 4.75 s. Furthermore, at 573 K, the crystallization of the surface of the deposited strand started at 12 s. The temperature range to keep the deposited strand in the amorphous

state was relatively narrow (15 K), and T_m approached 623 K. The mechanical properties of PEEK should be sufficiently high at these temperatures to permit stratification. At temperatures below 558 K, the kinetics of crystallization were fast, leading to crystallization before the deposition of a second strand, oppositely to the case at temperatures higher than 573 K. Consequently, the proposed environment temperature to promote interdiffusion was between 558 and 573 K. Furthermore, the temperature selection also depended on the printing velocity and strategy; specifically, our model allowed the optimization of the time between two depositions to keep the material in its amorphous phase.

Our results on crystallization by numerical simulation are in line with the experimental studies carried out by Kumar and Anderson [65] concerning the shape of the crystalline phases and Tardif et al. [59] concerning the peak of Nakamura kinetics for PEEK.

The observations by Tseng et al. [14] confirm our numerical results concerning the crystallization of PEEK. The fast crystallization of PEEK led to the highest degree of crystallinity of the final part, which was almost independent of the cooling rate of
the printed part. Although the printed parts were subjected to
heat treatment, there was no improvement in the mechanical
properties of the material. We simply suppose the influence of
the shrinkage of the final part. To further optimize the mechan-
ical properties, the experimental temperature for the liquefier
was 643–663 K and for the substrate up to 553 K, in line with our
numerical measurements.

Rinaldi et al. [66] printed PEEK samples at different cooling rates, reporting the mechanical properties of the printed parts under different conditions using tensile tests. The thermal transition investigated by DSC and the microstructure by X-ray diffraction and optical microscopy did not highlight any difference in terms of thermal transition, such as T<sub>g</sub>, T<sub>m</sub>, and final degree of crystallinity. The experimental results also confirmed that the degree of crystallinity in the printed parts was independent of the cooling rate of PEEK.

**6** | Conclusions

In this study, material deposition, heat transfer, and crystal-
lization in the MEX process were modeled using numerical
simulation. Furthermore, we determined the isothermal and
non-isothermal crystallization of PEEK by DSC experiments on
cooling for a large temperature range from the melting tempera-
ture up to the glass transition.

The half-time crystallization of PEEK 450G was determined experimentally by DSC for temperatures between 583 and 595 K. The experimental results were extrapolated using the Hoffman–Lauritzen equation to determine the half-time crystallization from T<sub>m</sub> (623 K) to T<sub>g</sub> (433 K). The comparison of the results for the Avrami kinetics coefficient by the predictive Hoffman–Lauritzen equation and Avrami equation shows very good agreement. Our results show that the crystallization rate of PEEK 450G is maximum at 512 K.

Finally, Avrami kinetics was applied to the Nakamura equation to determine the non-isothermal crystallization of PEEK 450G. By adding a convection term to the Nakamura non-isothermal crystallization, we determined the evolution of the crystallization field during the cooling of the semicrystalline polymers while considering the fluid flow. A modified Nakamura crystallization equation with a convection term was implemented in our model. Our outcome on the crystallization of PEEK is in line with several experimental studies reported in the literature.

The MEX process of PEEK with non-isothermal crystalliza-
tion was simulated to determine the kinetics of crystallization
during the material deposition.

Modeling the crystallization kinetics and heat transfer in the MEX process is essential for optimizing and controlling the interdiffusion of the polymer chains of the deposited strands. The model proposed for real material deposition in a two-phase system is a new approach for determining non-isothermal crystallization while considering fluid flow.

Providing a model enables the selection of the best printing pa-
rameters, such as temperature of the environment and substrate,

to optimize the interdiffusion of the molecular chains and conse-
quently improve the mechanical properties. However, the influ-
ence of the deposition of further layers remains to be included. The
deposition of a second layer increases the temperature of the first
deposited strand but is not sufficient to melt the crystalline phase.

The deposition strategy of layers on each other influences the interdiffusion of two adjacent layers, which is favored by the second layer deposition occurring before the cooling and crystallization of the first layer. During the selection of the printing parameters, raster strategy deposition, environment temperature, and relaxation time should be considered. It should be noted that the proposed approach in this study can also be applied to other polymer processes.

Acknowledgments

This study received financial support from the University of Toulouse, France, and Occitanie Region of France through a Ph.D. fellowship. Open access publishing facilitated by Tampereen yliopisto ja Tampereen ammattikorkeakoulu, as part of the Wiley - FinELib agreement.

Funding

This work was supported by Institut Clément Ader, Université de
Toulouse.

**Conflicts of Interest**

The authors declare no conflicts of interest.

Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

References

1. B. Hu, X. Duan, Z. Xing, et al., “Improved Design of Fused Deposition Modeling Equipment for 3D Printing of High-Performance PEEK Parts,” *Mechanics of Materials* 137 (2019): 103139, https://doi.org/10.1016/j.mechmat.2019.103139.

2. M. F. Arif, S. Kumar, K. M. Varadarajan, and W. J. Cantwell, “Performance of Biocompatible PEEK Processed by Fused Deposition Additive Manufacturing,” *Materials and Design* 146 (2018): 249–259, https://doi.org/10.1016/j.matdes.2018.03.015.

3. M. Doumeng, F. Berthet, K. Delbé, O. Marsan, J. Denape, and F. Chabert, “Effect of Size, Concentration, and Nature of Fillers on Crystallinity, Thermal, and Mechanical Properties of Polyetheretherketone Composites,” *Journal of Applied Polymer Science* 139, no. 5 (2022): e51574, https://doi.org/10.1002/app.51574.

4. R. Jones, P. Haufe, E. Sells, et al., “RepRap – The Replicating Rapid Prototyper,” *Robotica* 29 (2018): 177–191, https://doi.org/10.1017/S026357471000069X.

5. S. Upcraft and R. Fletcher, “The Rapid Prototyping Technologies,” *Assembly Automation* 23 (2007): 318–330.

6. K. M. Rahman, T. Letcher, and R. Reese, “Mechanical Properties of Additively Manufactured PEEK Components Using Fused Filament Fabrication,” in *Proceedings of the Proceedings of the ASME 2015 International Mechanical Engineering Congress and Exposition* (2015), V02AT02A009.

7. P. Honigmann, N. Sharma, B. Okolo, U. Popp, B. Msallem, and F. M. Thieringer, “Patient-Specific Surgical Implants Made of 3D Printed
PEEK: Material, Technology, and Scope of Surgical Application," *BioMed Research International* 2018 (2018): 1–8.

8. F. M. Thieringer, N. Sharma, A. Mootien, R. Schumacher, and P. Honigmann, "Patient Specific Implants From a 3D Printer – An Innovative Manufacturing Process for Custom PEEK Implants in Cranio-Maxillofacial Surgery," in *Industrializing Additive Manufacturing – Proceedings of Additive Manufacturing in Products and Applications-AMPA2017*, ed. M. Meboldt and C. Klahn (Springer International Publishing, 2018), 308–315.

9. X. Han, N. Sharma, Z. Xu, L. Scheideler, and J. Geis-gerstorfer, "An In Vitro Study of Osteoblast Response on Fused-Filament Fabrication 3D Printed PEEK for Dental and Cranio-Maxillofacial Implants," *Journal of Clinical Medicine* 8 (2019): 771, https://doi.org/10.3390/jcm8060771.

10. C. Basgul, T. Yu, D. W. MacDonald, R. Siskey, M. Marcolongo, and S. M. Kurtz, "Does Annealing Improve the Interlayer Adhesion and Structural Integrity of FFF 3D Printed PEEK Lumbar Spinal Cages?," *Journal of the Mechanical Behavior of Biomedical Materials* 102 (2020): 103455, https://doi.org/10.1016/j.jmbbm.2019.103455.

11. Y. Zhao, K. Zhao, Y. Li, and F. Chen, "Mechanical Characterization of Biocompatible PEEK by FDM," *Journal of Manufacturing Processes* 56 (2020): 28–42, https://doi.org/10.1016/j.jmapro.2020.04.063.

12. Y. Wang, W. Müller, A. Rumjahn, F. Schmidt, and D. Schwitalla, "Mechanical Properties of Fused Filament Fabricated PEEK for Biomedical Applications Depending on Additive Manufacturing Parameters," *Journal of the Mechanical Behavior of Biomedical Materials* 115 (2021): 104250, https://doi.org/10.1016/j.jmbbm.2020.104250.

13. F. Abdullah, K. Okuyama, A. Morimitsu, and N. Yamagata, "Effects of Thermal Cycle and Ultraviolet Radiation on 3D Printed Carbon Fiber Polyether Ether," *Aerospace* 7 (2020): 95, https://doi.org/10.3390/aerospace7070095.

14. J. W. Tseng, Y. Liu, K. Yen, et al., "Screw Extrusion-Based Additive Manufacturing of PEEK," *Materials & Design* 140 (2018): 209–221, https://doi.org/10.1016/j.matdes.2017.11.032.

15. M. F. Talbott, G. S. Springer, and L. A. Berglund, "The Effects of Crystallinity on the Mechanical Properties of PEEK Polymer and Graphite Fiber Reinforced PEEK," *Journal of Composite Materials* 21 (1987): 1056–1081, https://doi.org/10.1177/002199838702101104.

16. L. Martineau, F. Chabert, B. Boniface, and G. Bernhart, "Effect of Interfacial Crystalline Growth on Autohesion of PEEK," *International Journal of Adhesion and Adhesives* 89 (2019): 82–87, https://doi.org/10.1016/j.ijadhadh.2018.11.013.

17. B.-R. Cho and J. Kardos, "Consolidation and Self-Bonding in Poly(Ether Ether Ketone) (PEEK)," *Journal of Applied Polymer Science* 56 (1995): 1435–1454.

18. C. Yang, X. Tian, D. Li, Y. Cao, F. Zhao, and C. Shi, "Influence of Thermal Processing Conditions in 3D Printing on the Crystallinity and Mechanical Properties of PEEK Material," *Journal of Materials Processing Technology* 248 (2017): 1–7, https://doi.org/10.1016/j.jmatprotec.2017.04.027.

19. W. Z. Wu, P. Geng, J. Zhao, Y. Zhang, D. W. Rosen, and H. B. Zhang, "Manufacture and Thermal Deformation Analysis of Semicrystalline Polymer Polyether Ether Ketone by 3D Printing," *Materials Research Innovations* 18 (2014): S5-12–S5-16, https://doi.org/10.1179/1432891714Z.000000000898.

20. P. Wang, B. Zou, H. Xiao, S. Ding, and C. Huang, "Effects of Printing Parameters of Fused Deposition Modeling on Mechanical Properties, Surface Quality, and Microstructure of PEEK," *Journal of Materials Processing Technology* 271 (2019): 62–74, https://doi.org/10.1016/j.jmatprotec.2019.03.016.

21. S. Ding, B. Zou, P. Wang, and H. Ding, "Effects of Nozzle Temperature and Building Orientation on Mechanical Properties and

Microstructure of PEEK and PEI Printed by 3D-FDM," *Polymer Testing* 78 (2019): 105948, https://doi.org/10.1016/j.polymertesting.2019.105948.

22. M. Vaezi and S. Yang, "Extrusion-Based Additive Manufacturing of PEEK for Biomedical Applications," *Virtual and Physical Prototyping* 10 (2015): 123–135, https://doi.org/10.1080/17452759.2015.1097053.

23. M. Day, D. Sally, and D. M. Wiles, "Thermal Degradation of Poly (Aryl-Ether-Ether-Ketone): Experimental Evaluation of Crosslinking Reactions," *Journal of Applied Polymer Science* 40 (1990): 1615–1625, https://doi.org/10.1002/app.1990.070400917.

24. S. Bakrani Balani, H. Mokhtarian, E. Coatanéa, F. Chabert, V. Nassiet, and A. Cantarel, "Integrated Modeling of Heat Transfer, Shear Rate, and Viscosity for Simulation-Based Characterization of Polymer Coalescence During Material Extrusion," *Journal of Manufacturing Processes* 15 (2023): 108–114, https://doi.org/10.1016/j.jmapro.2023.02.021.

25. J. E. Seppala and K. D. Migler, "Infrared Thermography of Welding Zones Produced by Polymer Extrusion Additive Manufacturing," *Additive Manufacturing* 12 (2016): 71–76, https://doi.org/10.1016/j.addma.2016.06.007.

26. A. Costanzo, U. Croce, R. Spotorno, S. E. Fenni, and D. Cavallo, "Fused Deposition Modeling of Polyamides: Crystallization and Weld Formation," *Polymers* 12 (2020): 2980, https://doi.org/10.3390/polym12122980.

27. S. Bakrani Balani, F. Chabert, V. Nassiet, and A. Cantarel, "Influence of Printing Parameters on the Stability of Deposited Beads in Fused Filament Fabrication of Poly(Lactic) Acid," *Additive Manufacturing* 25 (2019): 112–121, https://doi.org/10.1016/j.addma.2018.10.012.

28. H. R. Vanaei, S. Fernandes, C. Fernando, M. Duarte, J. António, and M. Deligant, "Experimental Study of PLA Thermal Behavior During Fused Filament Fabrication," *Applied Polymer* 138, no. 4 (2020): 1–7, https://doi.org/10.1002/app.49747.

29. A. Bonito, M. Picasso, and M. Laso, "Numerical Simulation of 3D Viscoelastic Flows With Free Surfaces," *Journal of Computational Physics* 215 (2006): 691–716, https://doi.org/10.1016/j.jcp.2005.11.013.

30. T. Köpplmayr and E. Mayrhofer, "Numerical Simulation of Viscoelastic Layer Rearrangement in Polymer Melts Using OpenFOAM," *AIP Conference Proceedings* 1664 (2015): 1–5, https://doi.org/10.1063/1.4918407.

31. R. Comminal, J. Spangenberg, and J. H. Hattel, *Numerical Simulations of Viscoelastic Flows With Free Surfaces* (Congrès Français de Mécanique, 2013), 1–6.

32. D. Jun, W. Zhengying, W. Xin, W. Jijie, and C. Zhen, "An Improved Fused Deposition Modeling Process for Forming Large-Size Thin-Walled Parts," *Journal of Materials Processing Technology* 234 (2016): 332–341, https://doi.org/10.1016/j.jmatprotec.2016.04.005.

33. S. Dabiri, S. Schmid, and G. Tryggvason, "Fully Resolved Numerical Simulations of Fused Deposition Modeling," in *Proceedings of the Proceedings of the ASME 2014 International Manufacturing Science and Engineering Conference* (2014).

34. R. Comminal, *Numerical Simulation of Viscoelastic Free Surface Flows Using a Streamfunction/Log Conformation Formulation and the Volume of Fluid Method* (s, 2015).

35. R. Comminal, M. P. Serdeczny, D. B. Pedersen, and J. Spangenberg, "Numerical Modeling of the Strand Deposition Flow in Extrusion-Based Additive Manufacturing," *Additive Manufacturing* 20 (2018): 68–76, https://doi.org/10.1016/j.addma.2017.12.013.

36. R. Comminal, M. P. Serdeczny, D. B. Pedersen, and J. Spangenberg, "Motion Planning and Numerical Simulation of Material Deposition at
Corners in Extrusion Additive Manufacturing,” *Additive Manufacturing* 29 (2019): 100753.

37. A. D. Amico and A. M. Peterson, “An Adaptable FEA Simulation of Material Extrusion Additive Manufacturing Heat Transfer in 3D,” *Additive Manufacturing* 21 (2018): 422–430, https://doi.org/10.1016/j.addma.2018.02.021.

38. R. Jerez-Mesa, J. A. Travieso-Rodriguez, G. Gomez-Gras, X. Corbella, and R. Busqué, “Finite Element Analysis of the Thermal Behavior of a RepRap 3D Printer Liquefier,” *Mechatronics* 36 (2016): 119–126, https://doi.org/10.1016/j.mechatronics.2016.04.007.

39. S. Bakrani Balani, H. Mokhtarian, T. Salmi, and E. Coatanéa, “An Investigation of the Influence of Viscosity and Printing Parameters on the Extrudate Geometry in the Material Extrusion Process,” *Polymers* 15 (2023): 2202, https://doi.org/10.3390/polym15092202.

40. S. Bakrani Balani, A. Cantarel, F. Chabert, and V. Nassiet, “Influence of Parameters Controlling the Extrusion Step in Fused Filament Fabrication (FFF) Process Applied to Polymers Using Numerical Simulation,” *AIP Conference Proceedings* 1960 (2018): 140003.

41. F. Liravi, R. Darleux, and E. Toyserkani, “Additive Manufacturing of 3D Structures With Non-Newtonian Highly Viscous Fluids: Finite Element Modeling and Experimental Validation,” *Additive Manufacturing* 13 (2016): 113–123, https://doi.org/10.1016/j.addma.2016.10.008.

42. X. Li, Q. Wei, J. Li, et al., “Numerical Simulation on Crystallization-Induced Warpage of Injection-Molded PP/EPDM Part,” *Journal of Polymer Research* 26 (2019): 228.

43. B. Anass, M. Boutaous, R. El Otmani, et al., “Simulation of Crystallization Evolution of Polyoxymethylene During Microinjection Molding Cycle,” *Polymers for Advanced Technologies* 31, no. 4 (2019): 1–15, https://doi.org/10.1002/pat.4819.

44. R. Spina, M. Spekowius, and C. Hopmann. “Multiphysics Simulation of Thermoplastic Polymer Crystallization,” *Materials & Design* 95 (2016): 455–469, https://doi.org/10.1016/j.matdes.2016.01.123.

45. C. McIlroy and R. S. Graham, “Modelling Flow-Enhanced Crystallisation During Fused Filament Fabrication of Semi-Crystalline Polymer Melts,” *Additive Manufacturing* 24 (2018): 323–340, https://doi.org/10.1016/j.addma.2018.10.018.

46. Victrex, “Victrex PEEK 450g” (2014), 1–2.

47. N. Shibuya and R. S. Porter, “A Kinetic Study of PEEK Sulfonation in Concentrated Sulfuric Acid by Ultraviolet-Visible Spectroscopy,” *Polymer* 35 (1994): 3237–3242.

48. E. W. Lemmon and R. T. Jacobsen, “Viscosity and Thermal Conductivity Equations for Nitrogen, Oxygen, Argon, and Air,” *International Journal of Thermophysics* 25 (2004): 21–69.

49. N. Bond, “The Viscosity of Air,” *Proceedings of the Proceedings of the Physical Society* 49 (1926): 205.

50. S. Bakrani Balani, H. Mokhtarian, T. Salmi, and E. Coatanéa, “Numerical Simulation of the Shear Rate in the Fused Filament Fabrication Process of Poly-Ether-Ether-Ketone (PEEK),” *Progress in Additive Manufacturing* 10 (2024): 5057–5077, https://doi.org/10.1007/s40964-024-00888-3.

51. C. Bas, A. C. Grillet, F. Thimon, and N. D. Albérola, “Crystallization Kinetics of Poly(Aryl Ether Ether Ketone): Time-Temperature-Transformation and Continuous-Cooling-Transformation Diagrams,” *European Polymer Journal* 31 (1995): 911–921, https://doi.org/10.1016/0014-3057(95)00060-7.

52. P. Cebe and S. Hongt, “Crystallization Behaviour of Poly(Ether-Ether-Ketone),” *Polymer* 27 (1986): 1183–1192.

53. S. Tan, A. Su, J. Luo, and E. Zhou, “Crystallization Kinetics of Poly(Ether Ether Ketone) (PEEK) From Its Metastable Melt,” *Polymer* 40 (1999): 1223–1231.

54. J. Huang, “Dispersion, Crystallization Kinetics, and Parameters of Hoffman-Lauritzen Theory of Polypropylene and Nanoscale Calcium Carbonate Composite,” *Polymer Engineering and Science* 49, no. 9 (2009): 1855–1864, https://doi.org/10.1002/pen.

55. S. Vyazovkin, J. Stone, and N. Sbirrazzuoli, “Hoffman-Lauritzen Parameters for Non-Isothermal Crystallization of Poly (Ethylene Terephthalate) and Poly (Ethylene Oxide) Melts,” *Journal of Thermal Analysis and Calorimetry* 80 (2005): 177–180.

56. C. R. Choe and K. H. Lee, “Nonisothermal Crystallization Kinetics of Poly(Etheretherketone) (PEEK),” *Polymer Engineering and Science* (1989): 801–805.

57. C. N codeau, *Continuous Welding Modeling of Thermoplastic Matrix Composites* (Ecole Nationale Superieure d'Art et Metiers, 2005).

58. K. Nakamura, K. Katayama, and T. Watanabe, “Some Aspects of Nonisothermal Crystallization of Polymers. I. Relationship Between Crystallization Temperature, Crystallinity, and Cooling Conditions,” *Journal of Applied Polymer Science* 16 (1972): 1077–1091.

59. X. Tardif, B. Pignon, N. Boyard, et al., “Experimental Study of Crystallization of PolyEtherEtherKetone (PEEK) Over a Large Temperature Range Using a Nano-Calorimeter,” *Polymer Testing* 36 (2014): 10–19, https://doi.org/10.1016/j.polymertesting.2014.03.013.

60. E. Olsson, G. Kreiss, and S. Zahedi, “A Conservative Level Set Method for Two Phase Flow II,” *Journal of Computational Physics* 225 (2007): 785–807, https://doi.org/10.1016/j.jcp.2006.12.027.

61. D. Jacqmin, “Calculation of Two-Phase Navier–Stokes Flows Using Phase-Field Modeling,” *Journal of Computational Physics* 127 (1999): 96–127.

62. S. Osher and J. A. Sethian, “Fronts Propagating With Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formula tions,” *Journal of Computational Physics* 79 (1988): 12–49, https://doi.org/10.1016/0021-9991(88)90002-2.

63. L. D. Landau and E. M. Lifshitz, *Fluid Mechanics* (Pergamon Press, 1987).

64. S. Bakrani Balani, *Additive Manufacturing of the High-Performance Thermoplastics: Experimental Study and Numerical Simulation of the Fused Filament Fabrication* (University of Toulouse, 2019).

65. S. Kumar and D. P. Anderson, “Crystallization and Morphology of Poly (Aryl-Ether-Ether-Ketone),” *Polymer* 27 (1986): 329–336.

66. M. Rinaldi, F. Cecchini, L. Pigliaru, T. Ghidini, F. Lumaca, and F. Nanni, “Additive Manufacturing of Polyether Ether Ketone (PEEK) for Space Applications: A Nanosat Polymeric Structure,” *Polymers* 13, no. 11 (2021): 1–16, https://doi.org/10.3390/polym13010011.