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

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

# Integrated modeling of heat transfer, shear rate, and viscosity for simulation-based characterization of polymer coalescence during material extrusion

Shahriar Bakrani Balani<sup>a, *</sup>, Hossein Mokhtarian<sup>a</sup>, Eric Coatanéa<sup>a</sup>, France Chabert<sup>b</sup>, Valérie Nassiet<sup>b</sup>, Arthur Cantarel<sup>c</sup>

## ARTICLE INFO

### Keywords:

Coalescence

Material extrusion process

Fused filament fabrication

Finite element simulation

Polyether ether ketone

Polylactic acid

## ABSTRACT

The material extrusion process (MEX), also known as the fused filament fabrication process, has attracted attention in the manufacturing industry. A major obstacle to further application of the technology is the lack of mechanical strength due to the weak interlayer strength and poor coalescence between the adjacent beads. Understanding the effect of printing parameters on the coalescence of the adjacent beads is a step toward the improvement of the process. In this study, a novel two-phase flow numerical simulation approach coupled with heat transfer simulation has been applied to the high-viscosity polymers to determine the coalescence in the MEX process. The influence of printing temperature, substrate temperature, and the temperature of the printing chamber, as well as material deposition strategy (unidirectional and bidirectional) on the coalescence of the beads, has been investigated by numerical simulation and validated by experimental study. The modeling approach is applied to Glycerol, Polyether ether ketone (PEEK) and Polylactic acid (PLA). The results show that increasing temperature points (substrate temperature, chamber temperature, and printing temperature), increase the coalescence between the beads in the MEX process. The heat transfer model reveals that the cooling rate of the deposited bead in the MEX process is relatively high, and hence, the time window for reaching the coalescence between beads/layers is short. The heat transfer model also indicates that deposition of the further layers and beads does not influence the coalescence. The coalescence in the bidirectional deposition is higher compared to the unidirectional all conditions remaining similar. Unidirectional deposition leads to a uniform coalescence between the beads. However, the coalescence is not uniform for bidirectional deposition. The main novelty of this research is to simultaneously model heat transfer, shear rate and coalescence for numerical simulation to study the effect of printing parameters on the coalescence in the MEX process. Since the modeling of coalescence is time-consuming, two empirical equations based on obtained results have been proposed to predict the coalescence for PLA and PEEK separately.

### Nomenclature

$\vec{n}$ normal vector

$\dot{\gamma}$ shear rate [$s^{-1}$]

$\eta_0$ viscosity at terminal regime [Pa·s]

$\eta_{inf}$ viscosity at infinity shear rate [Pa·s]

$a$ dimensionless viscosity transition index [1]

$a_0$ initial radius of the filament [mm]

ABS Acrylonitrile butadiene styrene

$a_F$ final radius of the filament after coalescence [mm]

AM Additive Manufacturing

$a_T$ coefficient of Arrhenius law [1]

$C$ bonding length [mm]

$C_p$ specific heat capacity [$J·(kg·K)^{-1}$]

$d_n$ nozzle diameter [mm]

DAEs differential-algebraic system of equations

DSC differential scanning calorimeter

<footnote><sup>a</sup> Faculty of Engineering and Natural Sciences, Tampere University, Finland</footnote>

<footnote><sup>b</sup> LGP-ENIT-INPT, University of Toulouse, 47 Avenue d'Azereix, Tarbes, France</footnote>

<footnote><sup>c</sup> Institut Clément Ader (ICA), CNRS UMR 5312, University of Toulouse, IUT of Tarbes, UPS, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: shahriar.bakranibalani@tuni.fi (S. Bakrani Balani).</footnote>
<table><tr><td>E<sub>a</sub></td><td>activation energy [J·mol<sup>-1</sup>]</td></tr><tr><td>ε<sub>ls</sub></td><td>interfacial thickness parameter [mm]</td></tr><tr><td>F</td><td>other external forces [N]</td></tr><tr><td>FFR</td><td>filament feed rate [mm·s<sup>-1</sup>]</td></tr><tr><td>FFF</td><td>fused filament fabrication</td></tr><tr><td>F<sub>st</sub></td><td>surface tension force [kg·(m·s<sup>-2</sup>)]</td></tr><tr><td>g</td><td>gravity acceleration [m·s<sup>-2</sup>]</td></tr><tr><td>h</td><td>coefficient of convection [W·(m·K)<sup>-1</sup>]</td></tr><tr><td>H<sub>b</sub></td><td>height of the bead [mm]</td></tr><tr><td>H<sub>i</sub></td><td>layer height [mm]</td></tr><tr><td>k</td><td>coefficient of conduction [W·(m·K)<sup>-1</sup>]</td></tr><tr><td>L</td><td>length of the deposited bead [mm]</td></tr><tr><td>LS</td><td>level set</td></tr><tr><td>MEX</td><td>material extrusion process</td></tr><tr><td>n</td><td>pseudoplasticity index [1]</td></tr><tr><td>NBL</td><td>normalized bonding length [mm·mm<sup>-1</sup>]</td></tr><tr><td>ODEs</td><td>Ordinary Differential Equations</td></tr><tr><td>P<sub>a</sub></td><td>pressure applied in the liquefier [Pa]</td></tr><tr><td>P<sub>b</sub></td><td>pressure applied to the bead [Pa]</td></tr><tr><td>PEEK</td><td>Polyether ether ketone</td></tr><tr><td>PEKK</td><td>Polyether Ketone Ketone</td></tr><tr><td>PLA</td><td>Polylactic acid</td></tr><tr><td>P<sub>n</sub></td><td>pressure at the nozzle [Pa]</td></tr><tr><td>Q</td><td>heat source [W]</td></tr><tr><td>q<sub>Cond</sub></td><td>heat transfer by conduction [W·(m·K)<sup>-1</sup>]</td></tr><tr><td>Q<sub>e</sub></td><td>crystallization exothermic heat [W]</td></tr><tr><td>q<sub>Rad</sub></td><td>heat transfer by radiation [W·(m·K)<sup>-1</sup>]</td></tr><tr><td>R</td><td>universal gas constant [J·(mol·K)<sup>-1</sup>]</td></tr><tr><td>T</td><td>temperature [°C]</td></tr><tr><td>T<sub>0</sub></td><td>reference temperature for Arrhenius low [°C]</td></tr><tr><td>T<sub>amb</sub></td><td>temperature of ambient [°C]</td></tr><tr><td>T<sub>b</sub></td><td>temperature of the deposited bead [°C]</td></tr><tr><td>T<sub>c</sub></td><td>temperature of printing chamber [°C]</td></tr><tr><td>t<sub>d</sub></td><td>delay time [s]</td></tr><tr><td>T<sub>p</sub></td><td>printing temperature [°C]</td></tr><tr><td>TPF</td><td>two-phase flow</td></tr><tr><td>T<sub>s</sub></td><td>temperature of substrate [°C]</td></tr><tr><td>TS</td><td>travel speed [mm·s<sup>-1</sup>]</td></tr><tr><td>t<sub>t</sub></td><td>transition time (non-printing time) [s]</td></tr><tr><td>u</td><td>velocity field [mm·s<sup>-1</sup>]</td></tr><tr><td>W<sub>b</sub></td><td>width of the bead [mm]</td></tr><tr><td>γ</td><td>re-initialization parameter [m·s<sup>-1</sup>]</td></tr><tr><td>ε</td><td>coefficient of surface emissivity [1]</td></tr><tr><td>η</td><td>dynamic viscosity [Pa·s]</td></tr><tr><td>μ</td><td>surface tension [N·m<sup>-1</sup>]</td></tr><tr><td>ρ</td><td>density [kg·m<sup>-3</sup>]</td></tr><tr><td>σ</td><td>Stefan-Boltzmann constant [Kg·s<sup>-3</sup>·K<sup>-4</sup>]</td></tr><tr><td>ω</td><td>frequency [Hz]</td></tr><tr><td>θ</td><td>contact angle [Deg]</td></tr><tr><td>θ'</td><td>CONTACT angle revolution rate [Deg·s<sup>-1</sup>]</td></tr><tr><td>λ</td><td>relaxation time index [s]</td></tr><tr><td>φ</td><td>volume fraction [1]</td></tr></table>

# 1. Introduction

Additive Manufacturing (AM) of polymers is enabled mainly by material extrusion, powder sintering, and stereolithography technologies. Polymer additive manufacturing is based on the accumulation of successive layer depositions and creating the bonding through interlayer migration of the macromolecular chains through the interface. The material extrusion process (MEX) process is one of the wide varieties of polymeric AM processes and it is based on filament deposition. MEX is widely known as the fused filament fabrication (FFF) process. The principle of the FFF process is as follows: a solid filament enters the liquefier. The liquefier melts the filaments using a heating electric resistance [1]. The melted polymer is then deposited layer by layer,

according to a predefined trajectory. PLA (Polylactic acid) and, ABS (Acrylonitrile butadiene styrene) are among the most frequent polymers used in AM. Demanding applications in high-end fields such as aerospace and medical require high-performance polymers (with high mechanical properties). There is a limited number of high-performance polymers available commercially for the FFF process [2]. The mechanical and thermo-oxidation properties of Polyether ether ketone (PEEK) [3] make it a thermoplastic of significant interest for the AM processes such as the FFF and SLS (selective laser sintering) processes [4]. Since the introduction of PEEK to the market in 1980, abundant studies focused on material characterization and its applications in the high-performance domains [5–7]. However, it was only in 2015 that the first commercial FFF printer has been released to the market to print PEEK [8]. Nevertheless, contrary to the widely used polymers (e.g., PLA and ABS), printed parts with PEEK using the FFF process is not still well-mastered. This lack of specific knowledge of this material is an obstacle to reach the intended mechanical strength of the printed parts compared to the parts manufactured by traditional manufacturing processes. Wang et al. studied the influence of the nozzle diameter, nozzle temperature, and printing speed on the mechanical strength of PEEK parts printed by the FFF process [9]. They have shown that the nozzle diameter highly influences the mechanical strength of the PEEK parts.

Another process property that influences the mechanical strength of the part is the coalescence between the adjacent beads and layers [10]. Coalescence is the mechanism of bonding between two fluid droplets or two filaments. Mastering the coalescence phenomenon is essential in the manufacturing processes [11] such as sinter-based processing [12], and rotational molding [13]. The coalescence between beads determines the porosity of the parts. The optimization of the FFF process for high-quality parts requires mastering the influence of printing parameters on coalescence [14,15]. The kinetics of coalescence is described by a number of analytical models. Frankel [16] proposed the initial analytical model, which Eshelby later updated [17]. To regulate and enhance the coalescence, it is required to comprehend how temperature affects the rheological properties of the polymers. Despite its simplicity, it sheds light on how material characteristics affect the coalescence rate.

Le Boterf et al. studies the effect of thermal condition on coalescence kinetics through analytical study and experimental study for polypropylene (PP) [18]. They have shown that surface tension and viscosity are the most important driving forces of the coalescence. Their studies show that the effect of viscosity is more significant compared to the surface tension. Controlling the process parameters and modification of the polymer characteristics are among the most effective approaches to facilitate the macromolecular chain diffusion and entanglement between adjacent deposited beads and finally enhance the bonding interface of the beads [19]. Abbot et al. [15] investigated the influence of the process parameters (layer height, extruder temperature, and printing speed) on the coalescence of printed ABS parts. According to their findings, the printing speed significantly affects the tensile strength, which in turn depends on the coalescence between the deposited beads. Xu et al. [20] suggest that increasing the amorphous phase of the polymer structure could improve the interlayer bonding of additively manufactured PEEK parts. Their approach was to increase the amorphous state of the PEEK structure by compounding amorphous PAEK with PEEK blends (aPAEK/PEEK blends). Le Poivre et al. studied the surface tension and coalescence of PEKK (Polyether Ketone Ketone) by experimental study in isothermal and non-isothermal conditions [21]. They have studied the coalescence by the conjunction of a 2D heat transfer model with an analytical model [21]. Charlon et al. [22] investigated the influence of nozzle and chamber temperature on the adhesion and porosity of the printed PP parts. They have experimentally proved that increasing the chamber temperature improve the coalescence of the adjacent beads and adhesion of the beads.

According to the literature review performed by Al Rashid et al. on numerical simulation techniques for the FFF process [23], the modeling of coalescence is still a task requiring significant research efforts. The
two-phase flow (TPF) modeling is a powerful approach to model the complex systems when there are two phases in the system. It has widely been used by Comminal et al. to model the deposition of the material [24] and extrusion process [25,26]. The numerical simulation of coalescence has been mostly applied to low viscosity materials [27,28] in isothermal conditions [29]. Yuan et al. used the volume of fluid approach to model the coalescence of air bubble droplets and their breakup in a low viscosity continuous phase [30]. A recent study on the numerical simulation of the sintering was carried out by Rando et al. [31] in the food industry (chocolate BG-811). They have used the Navier-Stokes equations and Lagrangian-Eulerian approach to model the material flow through the interface. This approach has been applied to the material with a viscosity up to 700 Pa·s. Even though TPF numerical simulation approach is suitable for modeling the coalescence, in order to apply it to the FFF process it is required to determine the heat transfer during the material deposition. Consequently, measuring of the heat transfer in the FFF process predominates to determine the coalescence. Determination of the heat transfer in the FFF process has been the subject of several studies in recent years [32,33]. These studies are mainly based on experimental thermography [32,33], analytical studies, and numerical simulation [34,35]. Although the determination of the heat transfer by experimental study provides relatively accurate data, it is limited due to its costs and difficulty in operation. Experimentally determining the heat transfer is suitable for a perpendicular plane, and not across a 3D object. Consequently, numerical simulation is an alternative to overcome the limitations of the experimental studies.

Vanaei et al. [34] studied the influence of different process parameters (printing speed, temperature of substrate, and ambient temperature) on the heat transfer of a PLA printed wall using the FFF process by numerical simulation. The effect of the path planning and infill percentage on the heat transfer of the deposited bead by PLA has been studied by Ramos et al. [36]. Trofimov et al. [37] studied the heat transfer in the FFF process through experimental study and numerical modeling. Their studies demonstrated that numerical modeling is a powerful approach to computing the complex properties of the FFF process such as heat transfer and induced-temperature distortion. Their thermal studies on 3D printed parts in PLA show that the cooling rate during the FFF process is relatively high.

In order to simulate the temperature variation in the FFF process, Zhang et al. developed a 3D transient-state model [38]. Owens et al. [39] modeled the thermal history of the MEX processes by developing a 2D numerical model. A large number of thermal studies on the FFF process demonstrate its important influence on the properties of the printed parts. Although the current state of the art provides particularly good insight into the thermal condition of the FFF process, it is still limited to modeling a wall geometry or a single bead. The geometry of the wall or single bead does not provide information about the effect of deposition of several layers and beads on the heat transfer in the FFF process. In addition to temperature variation during the deposition, the flow behaviour of the polymer in the liquefier and after exiting from the nozzle highly influences the quality of the deposited bead and coalescence. The viscosity of the polymer depends on the flow characteristic such as pressure and velocity and liquefier geometry. Jayaraghul et al. studied the flow characteristics and heat transfer of PEEK used in the FFF process by experimental study and numerical simulation [40]. Mishara et al. studied the melt flow behaviour of PLA and ABS using numerical simulation [41]. In this study, the relationship between printing quality and pressure, temperature, velocity and viscosity inside the nozzle has been investigated. Serdeczny et al. [42] model material deposition on the substrate. The same approach could be employed in order to determine the viscoelastic behaviour of the polymer melt on the substrate. In our previous studies, we have used TPF numerical simulation approach in order to model the material extrusion and material deposition on the deposition substrate [8,43]. This approach could be used to model the material deposition and viscoelastic properties during material deposition.

The literature review lacks pertinent studies on the coalescence of the material with complex rheological properties. The published articles don't address the numerical simulation of the filament coalescence occurring during the FFF process. Consequently, the creation of a numerical simulation that considers the influence of deposition of several beads and layers on the heat transfer is essential to study the process properties related to the heat transfer such as coalescence, thermal warpage, and crystallization.

To better summarize the parameters/variables involved in the FFF process and their relationships, a colored directed graph is developed in Fig. 1 [44,45]. This graph describes the behavioral process taking place in the FFF process. Note that this oriented graph is not unique and different similar graphs can be developed depending on the study focus, intended details, and adopted assumptions as well as intended objectives. The variables in the directed graph are categorized into four principal groups with specified colors. The independent variables (green) are input process parameters or independent design variables that can be freely set. The independent variables are not causally dependent on other variables. Layer height ($H_l$), nozzle travel speed (TS), filament feed rate (FFR), and different temperature points, namely, printing temperature ($T_p$), substrate temperature ($T_s$), and the chamber temperature ($T_c$) are the independent variables in the scope of this research paper. The exogenous variables (grey) are the variables that are outside the system of interest, kept fixed in the study, or imposed on the system. The material properties as well as the nozzle-related variables are considered as the exogenous variables in the oriented-graph. Note that, even though the material properties change based on the temperature, they have been considered fixed values in this study as a modeling assumption. Dependent variables (blue) are causally dependent on other variables including independent and exogenous variables. The dependent variables are indirectly controllable. As it is shown in the oriented-graph, the temperature of the deposited bead ($T_b$) is the central variable in the graph. It is influenced by the set temperature points, crystallization exothermic heat ($Q_e$), and bead geometry variables ($W_b$, $H_b$) that are themselves influenced by the input process parameters (TS, $H_l$, and FF). The temperature of the deposited bead ($T_b$) influences dynamic viscosity ($\eta$) and surface tension ($\mu$). The bonding length (C) depends on polymer dynamic viscosity ($\eta$), surface tension ($\mu$), and the pressure applied to the deposited bead ($P_b$). In the current study, coalescence is the performance variable (red node). The coalescence is determined as normalized bonding length (NBL), which is the ratio between bonding

<figure><img src="image_3.png" /><figcaption>Fig. 1. Causal graph representation of the parameters influencing the coalescence in the FFF process.</figcaption></figure>
length (C) and nozzle diameter (d_n).

The ultimate objective of the current research is to improve the mechanical properties of the printed part using the FFF process by increasing the bonding quality of the adjacent beads. Increasing the bonding length of the adjacent deposited beads results in decreasing porosity; hence improving the mechanical properties. This objective is pursued by studying the effect of printing temperature (T_p), substrate temperature (T_s), and the temperature of the chamber (T_c) on the coalescence phenomena via numerical simulation and experimental validation. The coalescence and heat transfer numerical models are applied to PLA and PEEK as two materials of interest. The main novelty of this research is to simultaneously consider heat transfer and coalescence numerical modeling to study the effect of process parameters on the coalescence of adjacent beads in the FFF process. The novel coalescence numerical model is developed by applying level set (LS) equations to model the high viscosity fluids (such as PEEK). To go beyond the basic geometry of a wall or a single bead, in this research, we have proposed a matrix (which consists of five layers, and 4 adjacent beads in each layer).

The remainder of the manuscript is organized as the following: Section 2 describes the overall methodology and workflow used in this research. Section 2 also presents the characterization of the material used in this study. Section 3 describes the experimental setups. The development of finite element simulation and the predictive model used in this study are explained in Section 4. The results are presented and

discussed in Section 5. Finally, conclusions and research perspectives are addressed in Section 6.

## 2. Materials and methods

### 2.1. Methodology and workflow

The current study follows a step-by-step workflow shown in Fig. 2. In the first step, the material characterization has been performed for the target polymers of the study. In the second step, the modeling of coalescence for Glycerol was initiated with an experimental study. The boundary conditions derived from this experimental study are used to develop the numerical model of the coalescence for Glycerol using finite element simulation. Glycerol has a constant viscosity and surface tension at ambient temperature.

For Glycerol, the surface tension and viscosity are independent of temperature. This independency makes the implementation of the numerical modeling and experimental effort easier. The coalescence of the glycerol is much quicker compared to the polymers with relatively high viscosity. Despite this difference, a validated coalescence model for Glycerol serves as a basis to further develop the coalescence model for the target polymers of the study. Once the Glycerol coalescence model is validated, the numerical simulation is applied to model the coalescence of PLA and PEEK, separately. Similarly, coalescence models for PLA and

<figure><img src="image_4.png" /><figcaption>Fig. 2. Study workflow.</figcaption></figure>
PEEK are developed with experimental study followed by the development of finite element simulations, and analytical studies in step 3 and step 4. The coalescence numerical model for PLA and PEEK (simulation 2 in Fig. 2) is based on the validated numerical model for Glycerol (simulation 1 in Fig. 2), experimental settings (experiment 3 in Fig. 2), and the input of material characterization conducted earlier (experiment 1 in Fig. 2). This model is also validated with experimental and analytical studies.

Until the end of step 4, the material characterization, simulations and validating experiments have been performed in the test environment. In the test environment, two filaments are put together under controlled conditions with a regulated temperature increase. Hence, the coalescence simulations, model the coalescence of two adjacent filaments as it happens experimentally in the test environment. In order to determine the coalescence in the FFF process for a given geometry, it is essential to simulate the temperature variation according to the printing parameters, process environment conditions, and geometry. In step 5, the heat transfer models are developed to determine the temperature variation and associated viscosity for PLA and PEEK, considering the FFF process parameters in a given geometry. The final coalescence model for PLA and PEEK in the FFF process is then developed in Step 6 by inserting the results of the viscosity from the heat transfer model and shear rate in the coalescence model developed previously. Finally, as a model validation, the NBL value determined by numerical simulation has been compared to the NBL value of the part manufactured by the FFF process.

The driving modeling methodology guiding this research is based on the numerical simulation using the finite element methods and experimental validation using laboratory equipment. The experimental settings and numerical simulations are explained in detail in Sections 3 and 4 respectively.

## 2.2. Material characterization

To cover a broad range of polymers, PLA and PEEK have been chosen. PLA is widely used in the FFF process with relatively low viscosity. In contrast, PEEK is a high-performance thermoplastic with a relatively high viscosity, high melting temperature, and high mechanical properties. PLA and PEEK are both semi-crystalline materials which makes this

study comparable for both polymers. PEEK is suitable for high-performance demanding applications such as biomedical [46] [47] and aerospace applications [48]. PEEK 450G polymer pellet manufactured by Victrex and PLA NaturePlast PLI 005 polymer pellet were used in this study. PEEK filaments were produced in-house using a single screw extruder. To eradicate the moisture, the PEEK and PLA pellets were dried respectively for 24 h at 140 °C, and 3 h at 60 °C in the oven and stored in the desiccator [49].

As it is perceived from the causal graph shown in Fig. 1, surface tension and viscosity are directly influencing the coalescence rate of polymers. Surface tension and viscosity are material properties that are determined through test procedures. Determining the transition temperatures of polymers is a pre-requisite step for the determination of viscosity. First, the thermal transitions have been experimentally determined for the PLA and PEEK using a Q200–TA differential scanning calorimetry (DSC) apparatus under nitrogen gas flow. The DSC curves in Fig. 3 represent that the glass transition of PLA occurred at 60 °C, and the melting is completed at 175 °C. The crystallization of PLA in the cooling phase starts at 130 °C and ends at 100 °C. On the other side, the glass transition of PEEK is occurred at 150 °C, and the material is fully melted at 360 °C. The crystallization of PEEK starts from 310 °C and ends at 270 °C. The coalescence takes place in the cooling phase from the melting until when the crystallization starts. Hence, the coalescence takes place in the cooling phase from the printing temperature until 130 °C for PLA, and until 310 °C for PEEK.

The time-temperature equivalence indicates that the shear rate and temperature have significant impacts on the fluids' viscosity [50]. This is represented by Arrhenius law. As it is represented in Eq. (1), the viscosity ($\eta$) of the fluids is the shear rate ($\dot{\gamma}$) and temperature ($T$) dependent. Increasing the temperature and shear rate decreases the viscosity of the shear-thinning materials. Arrhenius coefficient ($a_T$) in Eq. (1) is determined by Eq. (3).

$$ \eta(\dot{\gamma}, T) = a_T \eta(\dot{\gamma}, T_0) \qquad (1) $$

The Carreau-Yasuda model [50], shown in Eq. (2), represents that the viscosity ($\eta$) depends on the shear rate ($\dot{\gamma}$) for non-Newtonian fluids. In Eq. (2), n represents the pseudoplasticity index, $\eta_0$ and $\eta_{inf}$ represent the viscosity of the fluid at the terminal regimes (zero shear rate and

<figure><img src="image_5.png" /><figcaption>Fig. 3. The representation of the thermal transition of PLA and PEEK by DSC curves.</figcaption></figure>
infinite shear rate). λ stands for the relaxation time index. ˙γ represents the shear rate and a represents the transition parameter.

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

Arrhenius law allows determining the viscosity at all temperature ranges above the crystallization. It is determined in the cooling phase from the melted state for a given material based on a known viscosity at a specified temperature. The Arrhenius coefficient (a_T) is determined by Eq. (3), where E_a indicates the activation energy [51] [52] and R is the universal gas constant [51].

$$ a_T = \exp \left[ \frac{E_a}{R} \left( \frac{1}{T} - \frac{1}{T_0} \right) \right] \quad (3) $$

In Arrhenius's law, the master curve represents the viscosity curve at T_0. a_T is equal to one at T = T_0. For the shear-thinning fluids and above T_0, a_T is greater than one (a_T > 1). This means that viscosity and temperature are inversely related.

Note that the time-temperature equivalence cannot be used in proximity to thermal transitions (the glass transition, crystallization temperature, and melting temperature). Furthermore, the Arrhenius law is not valid when the polymer is subjected to chemical structure change due to degradation, or evaporation. Our assumption is that the coalescence takes place at the amorphous state of the polymer. Consequently, the variation of the viscosity from the melted state until the start of the crystallization has been considered for the coalescence study. Furthermore, it is assumed that after exiting the melted polymer from the nozzle and during the coalescence there is no shear rate applied to the polymer (˙γ ≈ 0), consequently, the viscosity is at the terminal regime (η_0).

The complex viscosity from low frequency (Below 0.1 Hz) until high frequency (100 Hz) has been determined for PLA and PEEK at four temperature points above the melting temperature, by using the parallel-plate setting of the ARES rheometer [8,53]. The test experiments have been carried out at temperature points of 175 °C, 185 °C, 195 °C, and 205 °C for PLA and 350 °C, 366 °C, 383 °C, and 400 °C for PEEK. In order to predict the viscosity at the terminal regime below the reference temperature until the crystallization temperature and above the reference temperature the Arrhenius law has been applied to PLA at T = 205 °C and PEEK at T = 383 °C [53]. The experimental values of the viscosity for PLA and PEEK at the terminal regime (η_0), as well as the theoretical value predicted by the Arrhenius law, are shown in Fig. 4. The experimentally determined values for viscosity are in accordance

<figure><img src="image_6.png" /><figcaption>Fig. 4. Viscosity vs. Temperature for PLA and PEEK, determined experimentally by parallel-plate rheometer and by Arrhenius law – The experimental values are marked on the graph.</figcaption></figure>

with the predicted values by Arrhenius law. The maximum error observed is 9 % occurred in 195 °C for PLA and 8 % at 365 °C for PEEK. The surface tension of the melted polymer was estimated 0.028 N·m⁻¹ and 0.015 N·m⁻¹ for PLA and PEEK respectively [8], using the Parachor approximation equation [8,54]. Even though the Parachor equation overestimates the surface tension values at the low temperature, it is more accurate at the high temperature [21].

### 3. Experimental set-up

A Glycerol droplet is deposited with a 1.06 mm diameter syringe of Digi-drop instrument, on another droplet as it is shown in Fig. 9. Two adjacent filaments (with a diameter of 1.65 mm) are utilized for the coalescence tests of PLA and PEEK. A 200 cm³ oven has been used for heating the filaments. The filaments were placed between two metallic blocks to ensure maintaining contact between the filaments. To determine the temperature and heat rate accurately during the experiment, a K-type thermocouple was located between the two adjacent filaments, and the data was recorded with a Graphtec-GL240 data logger. The recording frequency of the thermocouple is 100 Hz with a measurement accuracy of ±1.5 Kelvin. To ensure homogenous temperature distribution, an average slow heating rate (0.019 °C.min⁻¹) was imposed to the polymers. This allows achieving quasi-isothermal conditions during the experiment. Images were captured from the cross-sections of the filament during the experiment by a high-speed camera which was aligned with the axis of the filaments. A high-speed camera with a 35 mm diameter lens and 1:1.4 magnification scale, has been used to capture the evolution of the coalescence at an exposure time of 5000 ms, and a frame rate of 50 images per second. To ensure the repeatability of the experiments, the experiments were repeated three times for both PLA and PEEK. In order to determine the shape of the deposited bead and voids in the manufactured part by 3D printing an X-ray tomography (CT-Scan) apparatus (Viscom X8050) has been used. The printed part is mounted on a rotating stage and the intensity of the projection images has been recorded at 1°. The repetition of the experiments has been stopped after three times, since the results showed good accordance.

### 4. Finite element simulation development

Developing the numerical modeling to describe the coalescence phenomena in the FFF process is complex due to the interdependencies of the process parameters, material properties, thermal history, and deposition trajectories of the process. Therefore, the current research work conducts a series of numerical simulations to step-by-step model and validates the coalescence in the FFF process ( workflow shown in Fig. 2). The modeling workflow established in this research work consists of developing a validated numerical simulation for describing the coalescence for Glycerol and using it as a basis to numerically develop the coalescence model for PLA and PEEK. A separate finite element simulation is then developed for heat transfer modeling of the FFF process on a given geometry. The ultimate coalescence model for PLA and PEEK in the FFF process is then developed by inserting the results of the temperature and viscosity from the heat transfer model in the coalescence model. The boundary conditions and governing equations of finite element simulations are described in the sections below. The simulations are implemented in COMSOL Multiphysics 6.0 software. Table 1 summarizes the material properties used in the development of finite element simulations for this study.

#### 4.1. Coalescence numerical simulation

A 2D-axisymmetry numerical simulation for Glycerol composed of two phases (air and Glycerol) has been developed. The axisymmetric condition is applied since the geometry, boundary condition, and material are symmetric with respect to the nozzle axis. The 2D-axisymmetry is an approach to reduce simulation efforts from a 3D model to
<table><caption>**Table 1**
Material properties used in the numerical modeling.</caption><thead><tr><th rowspan="2">Parameters</th><th rowspan="2">Symbol</th><th rowspan="2">Unit</th><th colspan="4">Material</th></tr><tr><th>Air</th><th>Glycerol</th><th>PLA</th><th>PEEK</th></tr></thead><tbody><tr><td>Density</td><td>ρ</td><td>kg·m<sup>-3</sup></td><td>1.145</td><td>1260</td><td>1250</td><td>1300</td></tr><tr><td>Coefficient of conduction</td><td>k</td><td>W·(m·K)<sup>-1</sup></td><td>0.025</td><td>N/A</td><td>0.13</td><td>0.28</td></tr><tr><td>Specific heat capacity</td><td>C<sub>p</sub></td><td>J·(kg·K)<sup>-1</sup></td><td>1005</td><td>N/A</td><td>1800</td><td>1800</td></tr><tr><td>Surface tension</td><td>μ</td><td>N·m<sup>-1</sup></td><td>N/A</td><td>0.062</td><td>0.025</td><td>0.018</td></tr><tr><td>Activation energy</td><td>E<sub>a</sub></td><td>J·mol<sup>-1</sup></td><td>N/A</td><td>N/A</td><td>106</td><td>90</td></tr><tr><td>Viscosity</td><td>η</td><td>Pa·s</td><td>1.8 × 10<sup>-5</sup></td><td>1.42</td><td></td><td></td></tr></tbody></table>

a 2D geometry. The modeling is then continued by developing a 2D numerical model composed of air and polymer (PLA and PEEK). The geometry (two adjacent filaments with a given length) allows us to define a condition in 2D to reduce computational effort. Note that an axisymmetric condition could not be applied for modeling PLA and PEEK, since the geometry is not a spherical droplet, but two adjacent filaments. Fig. 5 represents the boundary conditions of modeling space for the developed coalescence numerical simulations.

The coalescence numerical simulations are developed using the two-phase flow (TPF) modeling approach with the level set (LS) eq. LS method coupled with Navier-Stokes equations is an approach to track the interface of immiscible fluids [55]. For each phase in the modeling space, the Navier-Stokes equation and continuity equation are solved to determine the velocity field according to the properties of the fluids. The continuity and Navier-Stokes equations are represented in Eqs. (4) and (5).

$$ \nabla \cdot \mathbf{u} = 0 \qquad (4) $$

$$ \rho \frac{\partial u}{\partial t} = \nabla \cdot [-pI + \eta(T)(\nabla u + (\nabla u)^T)] + \rho g + F_{st} + F \qquad (5) $$

where the flow properties such as velocity, and pressure are shown as u and P. The ρ and η represent the fluid density and viscosity respectively.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Boundary condition of modeling space for coalescence numerical simulation using finite element method (Top: PLA and PEEK; Bottom: Glycerol).</figcaption></figure>

In the equations above, F<sub>st</sub> describes the force resulting from the surface tensions. The other external forces are represented by F and g shows the gravity acceleration. The transport equation inserted to the Navier-Stokes equations is represented in Eq. (6).

$$ \frac{\partial \varphi}{\partial t} + \nabla \cdot (u \varphi) = \gamma \nabla \cdot \left( \varepsilon_{ls} \nabla \varphi - \varphi(1-\varphi) \frac{\nabla \varphi}{|\nabla \varphi|} \right) \qquad (6) $$

where ϕ is the volume fraction, t stands for time, γ and ε<sub>ls</sub> represent the re-initialization parameter and regulating the interfacial thickness parameter. In the modeling space (shown in Fig. 5), ϕ = 1 represents the polymer phase and ϕ = 0 is air. For more accuracy, the interface thickness has been refined by reducing the ε<sub>ls</sub>. Viscosity and density of the modeling space in TPF are determined as a function of volume fraction, represented in Eqs. (7) and (8).

$$ \eta = \varphi \eta_{\text{polymer}} + (1 - \varphi) \eta_{\text{air}} \qquad (7) $$

$$ \rho = \rho \rho_{\text{polymer}} + (1 - \rho) \rho_{\text{air}} \qquad (8) $$

The same approach and equations have been applied in order to model the material deposition on the substrate and consequently determine the shear rate during material deposition on the substrate. The details of the boundary conditions and simulation are previously presented in the previous publications [8,56]. Since the shear rate depends on the printing parameters such as FFR, TS and, d<sub>n</sub> for the following studies on the coalescence, these parameters are considered constant (20 mm·s<sup>-1</sup>, 15 mm·s<sup>-1</sup> and 0.3 mm respectively).

## 4.2. Heat transfer numerical simulation

As it is represented in the colored directed graph (Fig. 1), temperature modeling is important to determine the coalescence kinetics in the FFF process. Consequently, prior to modeling the coalescence, it is necessary to determine the influence of printing conditions on the temperature at the reference point between the two adjacent beads. Fig. 6 represents the geometry considered for the heat transfer model. The geometry consists of a matrix with 5 layers and 4 beads in each layer. A contact point between bead number 10 and 11 in the middle of the matrix has been selected as a reference point for recording the temperature variation. The reference point is selected in such a way to reduce the influence of the substrate temperature and thermal conditions of the environment on the temperature of the reference point. This reference point makes it possible to determine the influence of several layers and beads deposition on the heat transfer. The geometry allows us to study the effect of the material deposition strategy (unidirectional and bidirectional deposition). The model determines the heat released by the crystallization, the influence of the printing temperature (T<sub>p</sub>), substrate temperature (T<sub>s</sub>), and chamber temperature (T<sub>c</sub>) on the heat transfer of the reference point. Furthermore, the printing speed (TS) and cross-section of the beads are assumed to be constant and the influence of temperature on the material properties has been neglected.

As it is represented in Eq. (9) the time interval between the deposition of two adjacent beads for a given point (delay time (t<sub>d</sub>)) depends on the length of the bead and printing speed. In Eq. (9), the transition time (t<sub>t</sub>) is 0.25 s for the bidirectional case, which corresponds to the pause time at the end of the bead deposition, prior to moving to the next starting position. For the bidirectional deposition, since there is no need for nozzle displacement along the length of the bead, the delay time (t<sub>d</sub>) is 0.25 s. The transition time (t<sub>t</sub>) for the unidirectional case is 1.5 s, which corresponds to the traveling time required from the extreme points at L = 50 mm to L = 0 mm as shown in Fig. 6.

$$ \left\{ \begin{array}{l} t_d = t_t + \frac{L}{TS} \\ t_t = 1.5 \quad \text{Unidirectional} \\ t_t = 0.25 \quad \text{Bidirectional} \end{array} \right. \qquad (9) $$
<figure><img src="image_8.png" /><figcaption>Fig. 6. A schematic representation of a 5 × 4 matrix consist of 20 beads, considered as geometry in temperature determination finite element simulation.</figcaption></figure>

The cross-section is considered as an elongated oval with a height equal to the nozzle diameter (d_n = 0.5 mm). From the empirical point of view, the height of the layer is normally equivalent to the nozzle diameter. Since during printing the melted polymer is pushed through the nozzle, it is considered that the polymer is spread on the previous layer. This makes a broader contact between the beads in the previous layers (a line) compared to the contact between two adjacent beads. Consequently, the results represented here are valid for bonding of two adjacent beads in the same layer and not between two successive layers. Before polymer deposition, material properties of the matrix are considered as air, during and after deposition, material properties change from air to polymer. For the transition from the air properties to the polymer properties, an ordinary differential equation (ODE) has been defined. Heat transfer is determined by a heat transfer equation represented in Eq. (10) [57]. The heat transfer model considers the influence of conduction (Eq. (11)), convection (Eq. (12)), radiation (Eq. (13)), and the released heat due to the crystallization of semi-crystalline polymer.

$$ \rho C_p \frac{\partial T}{\partial t} + \rho C_p u \cdot \nabla T + \nabla \cdot q = Q \quad (10) $$

$$ q_{\text{Cond}} = -k \nabla T \quad (11) $$

$$ - \vec{n} \cdot q_{\text{Conv}} = h(T_{\text{amb}} - T) \quad (12) $$

$$ - \vec{n} \cdot q_{\text{Rad}} = \varepsilon \sigma (T_{\text{amb}}^4 - T^4) \quad (13) $$

where C_p is the heat capacity, k represents the coefficient of conduction, ρ is density and Q represents the heat source. The coefficient of convection (h) in the ambient environment is considered as 10 W·(m·K)^-1, and u is zero. In the equation above, ε is surface emissivity, T_amb is the room temperature, and σ is the Stefan-Boltzmann constant [57].

For the heat transfer modeling, Ordinary Differential Equations (ODEs) and differential-algebraic system of equations (DAEs) (dode COMSOL module) has been used to model mesh activation. Dode module in COMSOL Multiphysics 6.0 enables the application of ODEs in modeling the gradual deposition of material and changing properties of the bead from air to polymer, as well as heat transfer in solid for modeling the heat transfer. To reach the numerical stability and increase efficiency in computational time, simulations are carried out by simultaneously computing Comsol's Dode and heat transfer in solid modules (strong coupling). The heat transfer numerical simulation is solved using the MUMPS solver. For the coalescence study, a phase initialization and a time-dependent study has been defined for the study. To simulate the fully coupled coalescence simulation PARDISO solver has been used. The time step for the coalescence simulation is set to be 0.01 s for the lower viscosity when kinetics of coalescence is fast, and 0.05 s for the higher viscosities when the kinetics of coalescence is slow. Fig. 7 shows the mesh applied to the geometries and the mesh sensitivity graph in modeling bonding length. Free tetrahedral with minimum mesh size of 0.1 mm is applied to the matrix structure. To reach more accurate numerical computation, and avoid exponentially increasing the computation time and data, the time step is set to 0.05 s when simulating near the reference point, and 0.15 s for the rest of the geometry. Adaptive mesh is applied to the models, meaning that the mesh is finer near the edges where more accuracy is required compared to the center of the bead. Note that the mesh size represented in Fig. 7 refers to the mesh size on the beads and not the substrate mesh size.

The size of the mesh in the LS modeling must be adapted to the parameters of the LS equation (ε_1s and γ). Consequently, using coarse mesh size leads to numerical errors, and simulation could not be completed. Therefore, the sensibility analysis of the mesh on the results represented in Fig. 7 shows similar results for mesh size smaller than 9E-3 mm. Above 9E-3 mm, the simulation could not be solved due to the numerical instability.

### 4.3. Analytical study- coalescence

The experimental study of coalescence is time-consuming and costly. It requires precise determination of the temperature during the process. Hence, fully relying on an experimental approach to determine the coalescence for all printing conditions (process parameters, and materials) is not cost-efficient. Consequently, investigation of other methodologies such as numerical simulation and analytical study, combined with an experimental approach, is necessary to achieve an accurate coalescence model. An analytical model for the determination of coalescence is investigated here. There are numerous models that can be used to explain particle coalescence. These models were initially created for the selective laser sintering procedure, which involves powder sintering [58]. The sintering of particles follows three steps: The particles stick together; a neck forms, but the particles continue to exist as distinct individuals; and finally, the neck between the particles expands until the particles lose their individual identity. In this article, modified Frankel's
<figure><img src="image_9.png" /><figcaption>Fig. 7. Representation of the meshed boundaries (I: Heat transfer simulation applied on matrix; II: Coalescence simulation applied on two adjacent beads) and the mesh sensitivity of the coalescence results.</figcaption></figure>

model proposed by Pokluda is used [59], shown in Eq. (14).

$$ \theta' = \frac{\mu}{\eta a_0} \frac{2^a \sin\theta \cos\theta (2 - \cos\theta)^{1-a}}{(1 - \cos\theta)(1 + \cos\theta)^{1-a}} \quad (14) $$

where $\theta'$ is the derivate of the angle evolution over time, $\mu$ is the surface tension, $a_0$ is the initial radius of the deposited beads, and $\eta$ represents the viscosity, $\theta$ is the angle between the intersecting lines between the line passes through the center and neck of the coalescence and the horizontal line (Fig. 8). Fig. 8 schematically represents the coalescence phenomena of two adjacent beads.

## 5. Results and discussion

The coalescence of a low viscosity fluid (Glycerol), and two high

<figure><img src="image_10.png" /><figcaption>Fig. 8. The representation of the coalescence between two filaments from the cross-section view.</figcaption></figure>

viscosity thermoplastics (PLA and PEEK) have been investigated in this study. The influence of the filament deposition strategy (unidirectional and bidirectional), printing temperature (Tp), chamber temperature (Tc), and substrate temperature (Ts) are studied by numerical simulation. Two sections below discuss the results of developed numerical simulations. The results of coalescence models are discussed by defining a dimensionless number to omit the influence of filament (nozzle) diameter on the coalescence. Normalized bonding length (NBL) is defined as the ratio between bonding length (C) and the diameter of the deposited bead (2 × a₀). In the MEX process, layer height (Hᵢ) determines the diameter of the deposited bead, hence it is replaced in the denominator of the Eq. (15). Normalized bonding length (NBL) as a dimensionless number is shown in Eq. (15).

$$ NBL = \frac{C}{2 \times a_0} = \frac{C}{H_i} \quad (15) $$

### 5.1. Coalescence

Implementing the numerical simulation for the Glycerol coalescence requires an accurate measurement of the dimension of droplets at the very beginning (t = 0 s). Dimensions of droplet and nozzle geometry serve as the initial geometry for the numerical simulation. Therefore, the experimental study of Glycerol coalescence had been carried out prior to the numerical simulation. As Fig. 9 represents, the initial modeling geometry is composed of a circle and a quasi-semi-circle with a diameter of 1.85 mm. The boundary conditions and governing equations used in the implementation of this simulation are explained in the previous section. As it is expected, the results of numerical simulation indicate that the coalescence of the Glycerol takes place considerably faster when compared to the coalescence of the high viscosity materials such as thermoplastics. Fig. 9 shows the results of the numerical model with corresponding experimental images. A complete coalescence from getting into contact (t = 0 ms), necking (t = 20 ms to t = 80 ms), and separation of the droplet from the nozzle occurred at 200 ms. The coalescence of Glycerol is fast mainly due to its low viscosity and high surface tension. The comparison of the numerical simulation results with the corresponding experimental studies shows very good agreement for the entire coalescence evolution from the necking to the
<figure><img src="image_11.png" /><figcaption>Fig. 9. Comparison of the results of numerical simulation and the experimental study for coalescence modeling of Glycerol.</figcaption></figure>

separation. Hence, the numerical simulation for the low viscosity fluid is validated. The next step is to evaluate the validation of the TPF numerical model for the high viscosity polymers.

The coalescence of Glycerol takes place in isothermal condition and ambient temperature. In contrast, the coalescence of thermoplastics occurs in non-isothermal condition. Therefore, it is necessary to control, measure, and record the temperature variation during the test. Similar to the Glycerol case, precise measurement is essential to accurately model the coalescence for PLA and PEEK. The experiments were repeated three times for PLA and PEEK filaments in order to ensure the repeatability of the experiments. Figs. 10 and 11 show the images of the evolution of the coalescence for PLA and PEEK in the experiments. The results confirm that the coalescence starts above the melting temperature for both PLA and PEEK. Increasing temperature leads to an increase

in the bonding length (C) between the adjacent filaments. The kinetics of coalescence of thermoplastics is slower due to the high viscosity near the melting temperature. The temperature interval for the coalescence of PEEK (350 °C to 380 °C) is wider; starting when the melting of the crystalline phase is completed, and the polymer is amorphous. Comparing the NBL values of PLA and PEEK shows that the coalescence between PLA filaments occurs much faster than in PEEK filaments. It is mainly because PLA has a lower viscosity and higher surface tension compared to PEEK. Since the spread of the filaments on the substrate is very low, the influence of gravity on the shape of the filaments has been neglected compared to the influence of the surface tension and viscosity.

While the determination of the non-isothermal coalescence for the polymers is costly, time-consuming, and prone to methodologic errors, determination of the isothermal coalescence at high temperature and

<figure><img src="image_12.png" /><figcaption>Fig. 10. Comparison of coalescence value determined by experimental study and numerical simulation for PLA.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Comparison of coalescence value determined by experimental study and numerical simulation for PEEK.</figcaption></figure>

coalescence in the FFF process by the experimental study is difficult due to test requirements. Consequently, it is necessary to investigate other effective methodologies. The experimental results are then compared to the modified Frankel's equation (Eq. (14)) and numerical simulation. To compare the modified Frankel's equation and numerical simulation with experimental results, the recorded temperature during the experimental study is transformed to viscosity at zero shear rate and inserted as the input in these models.

The two coalescence numerical models for PLA and PEEK are developed based on the validated Glycerol model according to the viscosity, filament diameter, and temperature recorded in experimental studies. Figs. 10 and 11 show the NBL values determined by the experimental study, numerical simulation, and modified Frankel's model for PLA and PEEK respectively. The horizontal axis at the bottom of the diagram represents the temperature of adjacent filaments, and the horizontal axis at the top of the diagram represents the time associated with the exact recorded temperature. Reporting the variation of the temperature (two horizontal axes) is of importance since the increase of the temperature between two adjacent filaments is different from the oven temperature at the same moment. Considering two horizontal axes simultaneously reveals the rate of temperature increase in the contact area of the two adjacent filaments. Note that the coalescence images taken from the experimental study have been subjected to image processing for improving the quality of the visualization. Furthermore, the colors of the PLA and PEEK were different from each other at the melting state (PLA is white and PEEK is black). Therefore, the images related to PLA are not meant to be compared to the images related to PEEK, due to different thresholds applied in the image processing. The main purpose of these images is to illustrate the validation of the models against the experimental study and show the coalescence evolution for PLA and PEEK separately. The comparison of the NBL values from the experimental study, modified Frankel's predictive model and numerical simulation shows good agreements in the results. However, there are still several discrepancies at some points on the curve. There is an important deviation near the melting temperature for both PLA (from $t = 0$ s to $t = 145$ s, and from $T = 154.5$ °C to $T = 158$ °C) and PEEK (from $t = 350$ s to $t = 450$ s, and from $T = 347.7$ °C to $T = 350$ °C). The deviations are mainly due to the initial contacting line between two filaments which has not been taken into account for Frankel's model.

Although Frankel's model provides relatively good results for the coalescence of two adjacent beads, the model considers only the surface tension and viscosity as driving forces for coalescence. Frankel's model is thus not valid for the more complex systems where there are other parameters such as printing pressure, and nozzle diameter influencing the coalescence.

In addition, since PEEK turns into black color in the melted state, it causes difficulties for image processing and measurement of the NBL value and consequently influences the measured results (the high-speed camera captures black and white images). Furthermore, there is a possible systematic error in the temperature measurement using a thermocouple. The selection of the LS parameters in the numerical simulation, particularly the re-initialization parameter ($\gamma$) highly influences the output. The determination of the maximum velocity of the fluid in the system prior to simulation is difficult and often done by trial and error. In the numerical simulation, the thickness of the air-polymer interface must be as fine as possible. This is feasible by refining the mesh size and interfacial thickness ($\epsilon_{ls}$) in the LS equation.

## 5.2. Heat transfer

At this stage, the coalescence models for PLA and PEEK have been developed and validated. In the next step, the aim is to determine the viscosity variation during the deposition in the FFF process via heat transfer simulation. The heat transfer numerical model has been validated by analyzing the results of the heat transfer with experimental data available in the literature [33].

Fig. 12 represents the comparison of the results of numerical simulation and the results of an experimental study in the literature [33]. To be able to compare the results and hence validate the heat transfer simulation, a numerical simulation is conducted for a multi-layer wall manufactured by the MEX process. Note that the simulation is applied on ABS as material to be able to match the condition of existing experimental results in the literature. The agreement between the simulation and experimental results validates the heat transfer modeling approach conducted in this research.

Fig. 13 represents the influence of the beads and layers deposition on the temperature variation of the reference point (between beads 10 and 11, see Fig. 6) in the matrix for PLA at $T_p = 200$ °C. The circled numbers
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Validation of heat transfer model with experimental results derived from literature on a multi-layer wall geometry [33].</figcaption></figure>

represented in Fig. 13 indicate the bead number (see Fig. 6). Hence, the circled numbers facilitate the analysis of the deposition time and the temperature variation of the reference point due to the deposition of any

specific bead number. The matrices on the top-right corner in Fig. 13 show the heat dissipation after deposition of bead 11 and 15 as two example scenarios. Note that the simulations use the material properties of polymer for previously deposited bead and air properties for the beads that are not yet deposited. For instance, bead 11 has only conduction heat transfer with bead 10 (left-side bead) and bead 7 (previous layer bead). The heat dissipation is governed by conduction to the lower bead and the previously deposited adjacent bead. Heat transfer is governed by convection, radiation, and conduction. This can be seen also visually that the heat dissipation in Fig. 13 is not symmetrical in different directions.

The numbers on the temperature curve represent the number of beads in the matrix (see Fig. 6). The cooling rate of the reference point is fast which is undesired for the coalescence of two beads. The temperature decreases from printing temperature to substrate temperature in approximately five seconds. According to the printing temperature, the temperature window for the coalescence varies from 0.35 s to 1.5 s in case of a higher printing temperature. The temperature-time curve derived from the numerical simulation shows that the addition of more beads and layers slightly increases the temperature of the reference point. The increase of the temperature due to adding layers is more influential than adding a bead in the same layer. It is mainly because of the low thermal conductivity of the polymers that makes the influence of bead deposition on temperature variation negligible. Deposition of layers 14 and 15 which are in contact with beads 10 and 11 have the highest influence on the temperature variation of the reference point. However, increased temperature by layer deposition is not sufficient to affect the coalescence of the beads at the reference point. Crystallization of the semi-crystalline thermoplastics is an exothermal reaction that leads to heat release during crystallization at crystallization temperature (130 °C for PLA and 300 °C for PEEK). Fig. 13 highlights that, the heat released by crystallization leads to approximately 1 °C to 2 °C increases in the average temperature.

As expected, increasing the chamber temperature decreases the cooling rate which is a favorable condition for increasing the coalescence. At $T_c = 20$ °C the temperature of the reference point decrease to approximately 80 °C before deposition of the next bead (bead 12, see Fig. 6). However, in the case of $T_c = 60$ °C, the temperature increases to approximately 105 °C before deposition of the bead 12. This difference in temperature highlights the important influence of increasing printing

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Temperature variation for PLA at $T_p = 200$ °C: influence of bead deposition on coalescence at reference point (The numbers on the curve represent the beads' number in the matrix).</figcaption></figure>
chamber temperature ($T_c$) and temperature substrate ($T_s$) to increase the coalescence of the beads during the FFF process. At a fixed printing temperature of $T_p = 200$ °C and $T_s = 60$ °C, the NBL values are calculated 0.43 and 0.476 for $T_c = 20$ °C and $T_c = 60$ °C respectively. The increase in the NBL value suggests that increasing the chamber temperature ($T_c$) improves the coalescence between the beads in the FFF process.

The fast-cooling rate of the reference point reveals that the influence of printing strategy or delay time between deposition of the adjacent beads is negligible on the coalescence in a unidirectional deposition. The reference point reaches the crystallization temperature before the deposition of the successive beads. The heat provided by the deposition of the next beads is not sufficient to increase the temperature for influencing the coalescence. The delay time between two adjacent beads for any points along bead length is the same for unidirectional deposition (Eq. (9)). Therefore, the coalescence is uniform along the bead length in the unidirectional deposition. However, in the case of bidirectional bead deposition, delay time influences the coalescence. In another word, the coalescence or NBL is not uniform along the bead length, since the delay time between deposition of two adjacent beads in a bidirectional deposition is not the same (minimum in one extreme point, and maximum in the other extreme point of the bead).

Fig. 14 shows the results of the simulation for analyzing the influence of the beads and layers deposition on temperature variation of the reference point for PEEK at $T_p = 400$ °C. Similar to Fig. 13, the circled numbers represent the bead number in the matrix (see Fig. 6). As it was expected, similar conclusions are drawn from the results of both PLA and PEEK. The pattern in temperature variation in the matrix for PEEK is similar to the PLA, however, there are some differences between the curves. The cooling rate for PEEK is higher than PLA due to the higher conductivity of the PEEK compared to PLA. This leads to a lower temperature peak for PEEK during the deposition of bead 14 and bead 15 compared to PLA. The effect of exothermic heat released due to crystallization is negligible.

Two points near the extreme points of the matrix and reference points (between beads 10 and 11) are selected to investigate the influence of bidirectional deposition on the NBL value. In order to reduce the influence of the environment temperature on the selected points at the

extremities, the points were selected at a distance of 2 mm from the extreme points of the bead (Fig. 15:top). The obtained results of the viscosity and NBL values for the bidirectional deposition of PEEK are represented in Fig. 15. Note that viscosity represented in Fig. 15 is shear rate – temperature dependent. The shear rate of the polymer melt after exiting the nozzle is relatively high (in this specific case it is >200 s⁻¹) and since the thermoplastics are shear-thinning materials this highly reduces the viscosity of the deposited bead. The influence of the shear rate on the viscosity from zero to 0.25 s after the polymer exits from the nozzle is very important. The shear rate leads to a drastic drop in the viscosity immediately during deposition. In the case of PEEK, the combination of the shear rate and temperature results in reducing the viscosity below 100 Pa·s. However, the observations derived from simulation show that the influence of the shear rate on viscosity for PEEK is more important than PLA. This is mainly related to the

<figure><img src="image_16.png" /><figcaption>Fig. 14. Temperature variation for PEEK at $T_p = 400$ °C and $T_c = 158$ °C: influence of bead deposition on coalescence at the reference point (The numbers on the curve represent the same beads' number in the matrix shown in Fig. 13).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. The variation of the viscosity of three points between bead 10 and 11 for bidirectional deposition for $T_p = 500$ °C, $T_s = 150$ °C and $T_c = 20$ °C.</figcaption></figure>
viscoelastic behaviour of the PEEK compared to PLA and the parameters of the Carreau-Yasuda equation (represented in Fig. 4). In the zoomed curve represented in Fig. 15 there is a small jump (increase and then decrease) in the viscosity (from of 5 to 10 Pa·s) due to the influence of the extruder body on the deposited bead. By removing this pressure from the extruder polymer is returned to a release state. However, the influence of this variation is very slim and negligible.

Contrary to the unidirectional deposition, in the bidirectional deposition, the NBL value highly depends on the delay time. For the shorter delay time, the NBL value is higher than the longer delay time. In Fig. 15, the delay time for L = 48 mm is 1 s, while for L = 25 mm and L = 2 mm the delay time is 3.9 s and 6.95 s respectively. At the shorter delay time, the temperature difference between the printing temperature and the temperature of the previously deposited bead is lower compared to the longer delay time, which results in a slower cooling rate. This variation in the cooling rate results in viscosity variation which influence the coalescence in the bidirectional deposition. At L = 48 mm point, the maximum NBL is 0.55, while for the L = 2 mm the maximum NBL reaches 0.5. The obtained results for the coalescence in the bidirectional deposition recommend that to increase the coalescence between the beads, it is necessary to reduce the delay time during deposition (increasing the printing speed and reducing the length of the beads). Although theoretically coalescence takes place from the printing temperature (T<sub>p</sub>) until crystallization temperature, above the η = 5000 Pa·s, the kinetics of coalescence is very slow.

Fig. 16 shows the influence of the printing chamber temperature (T<sub>c</sub>) and substrate temperature (T<sub>s</sub>) on the NBL value of PEEK at T<sub>p</sub> = 500 °C. As expected, increasing the printing chamber temperature (T<sub>c</sub>) and substrate temperature (T<sub>s</sub>) highly increases the NBL value in the bidirectional deposition. The effect of printing chamber temperature (T<sub>c</sub>) on the NBL value is greater than the effect of substrate temperature (T<sub>s</sub>), due to the polymer's low conductivity (in this case, the coefficient of conduction is 0.13 W·(mK)<sup>-1</sup> for PEEK). Low thermal conductivity leads to low heat transfer from the substrate to the beads by heat conduction.

The linear trendline fitted on the NBL value of the PEEK in Fig. 16 for the bidirectional deposition case confirms that the variation of the NBL value in bidirectional is due to the difference in the cooling rate and the

<figure><img src="image_18.png" /><figcaption>Fig. 16. Influence of printing chamber temperature on the NBL value for the bidirectional deposition at T<sub>p</sub> = 500 °C.</figcaption></figure>

initial temperature before deposition of the adjacent bead at the selected points.

The results of simulations indicate that to increase the coalescence (NBL value), the printing temperature (T<sub>p</sub>) is among the most influential parameters. However, an excessive increase in printing temperature (T<sub>p</sub>) leads to polymer degradation. In contrast, increasing the printing temperature of chamber (T<sub>c</sub>) reduces the risk of polymer degradation, even though it may increase the printing and printer costs. The exposure time of the polymer to the high temperature must be selected according to the polymer thermal degradation. Consequently, the liquefier design can be optimized to reduce the risk of thermal degradation while increasing the printing temperature. In Fig. 16, the NBL values for T<sub>c</sub> = T<sub>s</sub> = 150 °C match perfectly with the NBL values of T<sub>c</sub> = 20 °C and T<sub>s</sub> = 250 °C. This highlights an alternative of increasing substrate temperature (T<sub>s</sub>) instead of having a printer with a regulated temperature of chamber. The alternative of increasing substrate temperature is more cost-effective and reduces the design complexity of the printer. However, the alternative of developing a machine with a closed chamber to regulate the temperature of chamber remains a possible approach to enhance the coalescence.

Fig. 17 depicts the NBL for the PLA and PEEK values at the reference point for unidirectional deposition according to different printing temperatures. Furthermore, the NBL of PLA is more sensitive to temperature increase (see curves in Fig. 17). For the temperature near to the melting temperature (For PLA: 180 °C and for PEEK: 350 °C) the coalescence is slower, while at the higher temperature the kinetics of coalescence increases.

The NBL value of PEEK at (T<sub>m</sub> = 500 °C, T<sub>c</sub> = 20 °C, and T<sub>s</sub> = 150 °C) in the reference point for unidirectional is 0.5 and for bidirectional deposition is 0.52. The effect of the gravity on the shape of the PLA beads is considerably high at T<sub>m</sub> = 220 °C. This is not observed for the PEEK at the higher temperature. This is mainly because of the surface tension of PLA which is higher than PEEK. The NBL values at the reference point for both unidirectional and bidirectional cases are the same because the delay time (t<sub>d</sub>) is longer than the time required for reaching the crystallization temperature. The NBL value is uniform for all the given points in the unidirectional deposition. Therefore, the NBL value for the unidirectional case can be indirectly determined by the simulation of bidirectional deposition. This is beneficial since the numerical simulations of the coalescence and heat transfer are time-consuming and computationally demanding.

Furthermore, since these simulations are time-consuming and computationally expensive, a simplified equation for the determination of the NBL value as a function of temperature is necessary. Such types of equations can be used for modeling the process on a bigger scale to create a meta-model for modeling the FFF process and process optimization. Consequently, Eqs. (16) and (17) show two polynomial curves fitted on the determined NBL value of PLA and PEEK (represented in Fig. 17) to predict the NBL value at the reference point according to the printing temperature.

$$NBL_{PLA} = (1*10^{-4}) * T_p^2 - 0.0366 * T_p + 3.3021 \quad (16)$$

$$NBL_{PEEK} = (1*10^{-5}) * T_p^2 - 0.0077 * T_p + 1.7057 \quad (17)$$

where T<sub>p</sub> is the printing temperature. R-squared (R<sup>2</sup>) for fitted results on NBL of PLA and PEEK is R<sup>2</sup><sub>PLA</sub> = 0.9946 and R<sup>2</sup><sub>PEEK</sub> = 0.9975, knowing that R<sup>2</sup> close to 1 shows the low dispersity of the fitted curve on the NBL results.

### 5.3. Experimental validation of coalescence model in FFF process

To experimentally assess the validity of the developed coalescence model, a large matrix geometry of 29 layers and 7 adjacent beads in each layer was printed. The sample is printed using a nozzle with a diameter of 0.3 mm. As printing parameters, layer height (H<sub>i</sub>), filament feed rate
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Coalescence of PLA and PEEK according to the printing temperature for unidirectional printing ($T_c$ = 20 °C and $T_s$ = 60 °C for PLA; $T_c$ = 20 °C and $T_s$ = 150 °C for PEEK).</figcaption></figure>

(FFR), and travel speed (TS) are set to 0.35 mm, 20 mm/s, and 30 mm/s, respectively. The printing temperature ($T_p$) and substrate temperature ($T_s$) are fixed at 200 °C and 60 °C, respectively. Fig. 18 shows an image of the cross-section of the printed matrix using an X-ray Tomography. This image allows the measurement of the bonding length (C) between two adjacent beads. The bonding length (C) of a total of N = 174 data points was measured using ImageJ software. It was observed for three points that the adjacent beads were not bonded together (C = 0). These points were discarded from the further analysis, even though their existence was not influencing the results due to the relatively high number of data points (N = 174). The normalized bonding length (NBL) is hence calculated according to the average actual bead height (0.28 mm) which is measured from the image. A quantile-quantile (Q-Q) plot is used for

evaluating the distribution of the measured bonding length (C) and hence calculated normalized bonding length (NBL). The Q-Q plot shown in Fig. 18 shows a linear relationship between measured quantiles (Z-score) with normal theoretical quantiles (Z-score) with $R^2$ = 0.992. The Q-Q plot confirms that the NBL follows a normal distribution with an average value (μ) of 0.52 and a standard deviation (σ) of 0.13 for NBL (N (0.49, 0.12)). A simulation is run with similar parameter settings as this experimental test. The simulation predicts NBL = 0.43 with an error of 11.51 % compared to the average experimental NBL = 0.49. The predicted value is the interval of (μ-σ, μ + σ) where it includes 64.36 % of the total measured data (112 data points).

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Experimental validation of coalescence model in FFF process and statistical indicators.</figcaption></figure>
## 6. Conclusion

In this study, the coalescence study of Glycerol, PLA, and PEEK has been performed by numerical simulation and experimental study. A novel TPF-LS numerical simulation approach has been applied to the high viscosity polymers. The numerical simulations have been validated with the experimental studies. In the current research, the TPF numerical simulation has been combined with the heat transfer and shear rate model to determine the coalescence of the two adjacent beads according to the printing parameters, for the FFF process.

The effect of the printing temperature, substrate temperature, and the temperature of the printing chamber, as well as the material deposition strategy (unidirectional and bidirectional), have been investigated. The results show that increasing temperature points, mentioned above, increase the coalescence between beads in the FFF process. Increasing the printing temperature is the most effective approach to increase the coalescence. However, increasing printing temperature increases the risk of polymer degradation. Therefore, the printing temperature setting is limited to the polymer thermal degradation (maximum 500 °C for PEEK and 220 °C for PLA). The increase in the temperature of printing chamber requires the printing in a regulated heating chamber which eventually increases the processing costs, the complexity of the printer, and printer cost. This makes the printing chamber temperature of the least favorable alternative to enhance the coalescence from the cost point of view.

The heat transfer numerical simulation shows that the cooling rate of the deposited bead in the FFF process is relatively fast, and hence, the time window for reaching the coalescence between beads/layers is short. The heat transfer simulation indicates that deposition of the further layers and beads does not influence the coalescence. The coalescence in the bidirectional deposition is higher compared to unidirectional in the same condition. Unidirectional deposition leads to a uniform coalescence between beads; however, the coalescence is not uniform for the bidirectional deposition. This highlights the conclusion that bidirectional deposition leads to better coalescence compared to unidirectionally and perimeter raster deposition strategies. In the case of perimeter raster deposition, the coalescence is lower compared to unidirectional and bidirectional deposition due to the longer delay time for perimeter raster deposition.

Although numerical simulation is more convenient than experimental study, it is still time-consuming and computationally demanding. For instance, it takes approximately 20 h for each heat transfer simulation and 4 h for each coalescence simulation (with a desktop computer Core i7-7700HQ CPU @ 2.80 Hz and 32 GB RAM). This highlights the need to develop a meta-model to predict the coalescence and optimize the FFF process.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This work was supported by the Finnish Institute of Technology (FITech) [grant number 310290011131], and Université Fédérale de Toulouse and Occitanie region of France through a Ph.D. fellowship APR2015.

## References

[1] Upcraft S, Fletcher R. The rapid prototyping technologies. Assem. Autom. 2007;23: 318–30.

[2] Golbang A, Harkin-Jones E, Wegrzyn M, Campbell G, Archer E, McIlhagger A. Production and characterization of PEEK/IF-WS2 nanocomposites for additive manufacturing: simultaneous improvement in processing characteristics and

material properties. Addit Manuf 2020;31:100920. https://doi.org/10.1016/j.addma.2019.100920.

[3] Kurtz S, editor. PEEK biomaterials handbook. Second: Elsevier; 2019.

[4] Schmidt M, Pohle D, Rechtenwald T. Selective laser sintering of PEEK. CIRP Ann - Manuf Technol 2007;56:205–8. https://doi.org/10.1016/j.cirp.2007.05.097.

[5] Zhang J, Tian W, Chen J, Yu J, Zhang J, Chen J. The application of polyetheretherketone (PEEK) implants in cranioplasty. Brain Res Bull 2019;153: 143–9. https://doi.org/10.1016/j.brainresbull.2019.08.010.

[6] Rinaldi M, Cecchini F, Pigliaru L, Ghidini T, Lumaca F, Nanni F. Additive manufacturing of polyether ether ketone (PEEK) for space applications: a nanosat polymeric structure. Polymers (Basel) 2021;13:11. https://doi.org/10.3390/polym13010011.

[7] Spece H, Yu T, Law AW, Marcolongo M, Kurtz SM. 3D printed porous PEEK created via fused filament fabrication for osteoconductive orthopaedic surfaces. J Mech Behav Biomed Mater 2020;109. https://doi.org/10.1016/j.jmbbm.2020.103850.

[8] Bakrani Balani S. Additive manufacturing of the high-performance thermoplastics: experimental study and numerical simulation of the fused filament fabrication. University of Toulouse; 2019.

[9] Wang Y, Müller WD, Rumjahn A, Schmidt F, Schwitalla AD. Mechanical properties of fused filament fabricated PEEK for biomedical applications depending on additive manufacturing parameters. J Mech Behav Biomed Mater 2021;115: 104250. https://doi.org/10.1016/j.jmbbm.2020.104250.

[10] Qi S, Gao X, Su Y, Dong X, Cavallo D, Wang D. Correlation between welding behavior and mechanical anisotropy of long chain polyamide 12 manufactured with fused filament fabrication. Polymer (Guildf) 2021;213:123318. https://doi.org/10.1016/j.polymer.2020.123318.

[11] Randall RG. Sintering theory and practice. 1996.

[12] Taki K, Tabata K, Kihara S-I, Ohshima M. Bubble coalescence in polymeric foaming process. In: Foam. 2004 4th Int. Conf. Thermoplast. Foam Process. Technol. - Tech. Pap; 2004. p. 1–39.

[13] Asgarpour M, Bakir F, Khelladi S, Khavandi A, Tcharkhtchi A. Characterization and modeling of sintering of polymer particles. J Appl Polym Sci 2011;119:2784–92. https://doi.org/10.1002/app.32924.

[14] Reffas SA, Elmeguenni M, Benguediab M. Analysis of void growth and coalescence in porous polymer materials. In: Coalescence in polymer materials. Eng. TechnolAppl. Sci. Res. 3; 2013. p. 452–60.

[15] Abbott AC, Tandon GP, Bradford RL, Koerner H, Baur JW. Process-structure-property effects on ABS bond strength in fused filament fabrication. Addit Manuf 2018;19:29–38. https://doi.org/10.1016/j.addma.2017.11.002.

[16] Frankel J. Viscous flow of crystalline bodies under the action of surface tension. J. Physics, USSR 1945;5:385.

[17] Eshelby JD. Discussion in paper by shaler AJ, seminar on kinetics of sintering. Met Trans 1949;185:806–7.

[18] Le Boterff J, Charlon S, Soulestin J. How does the temperature of polymer bead influence the kinetics of coalescence in additive manufacturing processes? J Appl Polym Sci 2022;139:1–11. https://doi.org/10.1002/app.51782.

[19] Gao X, Qi S, Kuang X, Su Y, Li J, Wang D. Fused filament fabrication of polymer materials: a review of interlayer bond. Addit Manuf 2021;37:101658.

[20] Xu Q, Xu W, Yang Y, Yin X, Zhou C, Han J, Li X, Shang Y, Zhang H. Enhanced interlayer strength in 3D printed poly (ether ether ketone) parts. Addit Manuf 2022;55:102852. https://doi.org/10.1016/j.addma.2022.102852.

[21] Lepoivre A, Levy A, Boyard N, Gaudefroy V, Sobotka V. Coalescence in fused filament fabrication process: thermo-dependent characterization of high-performance polymer properties. Polym Test 2021;98:107096. https://doi.org/10.1016/j.polymtesting.2021.107096.

[22] Charlon S, Le Boterff J, Soulestin J. Fused filament fabrication of polypropylene: influence of the bead temperature on adhesion and porosity. Addit Manuf 2021;38: 4–11. https://doi.org/10.1016/j.addma.2021.101838.

[23] Al Rashid A, Koç M. Fused filament fabrication process: a review of numerical simulation techniques. Polymers (Basel) 2021;13. https://doi.org/10.3390/polym13203534.

[24] Comminal R, Serdeczny MP, Pedersen DB, Spangenberg J. Numerical modeling of the strand deposition flow in extrusion-based additive manufacturing. Addit Manuf 2018;20:68–76. https://doi.org/10.1016/j.addma.2017.12.013.

[25] Comminal R, Serdeczny MP, Pedersen DB, Spangenberg J. Motion planning and numerical simulation of material deposition at corners in extrusion additive manufacturing. Addit Manuf 2019;29:100753.

[26] Serdeczny MP, Comminal R, Mollah T, Pedersen DB, Spangenberg J. Numerical modeling of the polymer flow through the hot-end in filament-based material extrusion additive manufacturing. Addit Manuf 2020;36:101454.

[27] Sellier M, Trelluyer E. In: Modeling the coalescence of sessile droplets; 2009. p. 1–14. https://doi.org/10.1063/1.3154552.

[28] Mohammadi MM, Shahhosseini S, Bayat M. Direct numerical simulation of water droplet coalescence in the oil. Int J Heat Fluid Flow 2012;36:58–71. https://doi.org/10.1016/j.ijheatfluidflow.2012.04.001.

[29] Xu Z. Numerical simulation of the coalescence of two bubbles in an ultrasound field. Ultrason Sonochem 2018;49:277–82. https://doi.org/10.1016/j.ulsonch.2018.08.014.

[30] Yuan S, Dabirian R, Shoham O, Mohan RS. Numerical simulation of liquid droplet coalescence and breakup. J Energy Resour Technol Trans ASME. 2020;142:1–11. https://doi.org/10.1115/1.4046603.

[31] Rando P, Ramaoli M. Numerical simulations of sintering coupled with heat transfer and application to 3D printing. Addit Manuf 2022;50:102567. https://doi.org/10.1016/j.addma.2021.102567.
[32] Ferraris E, Zhang J, Van Hooreweder B. Thermography based in-process monitoring of fused filament fabrication of polymeric parts. CIRP Ann 2019;68: 213–6. https://doi.org/10.1016/j.cirp.2019.04.123.

[33] Lepoivre A, Boyard N, Levy A, Sobotka V. Heat transfer and adhesion study for the FFF additive manufacturing process. Procedia Manuf 2020;47:948–55. https://doi.org/10.1016/j.promfg.2020.04.291.

[34] Vanaei HR, Khelladi S, Deligant M, Shirinbayan M, Tcharkhtchi A. Numerical prediction for temperature profile of parts manufactured using fused filament fabrication. J Manuf Process 2022;76:548–58. https://doi.org/10.1016/j.jmapro.2022.02.042.

[35] Vanaei HR, Shirinbayan M, Deligant M, Khelladi S, Tcharkhtchi A. In-process monitoring of temperature evolution during fused filament fabrication: a journey from numerical to experimental approaches. Thermo. 2021;1:332–60. https://doi.org/10.3390/thermo1030021.

[36] Ramos N, Mittermeier C, Kiendl J. Experimental and numerical investigations on heat transfer in fused filament fabrication 3D-printed specimens. Int J Adv Manuf Technol 2022;118:1367–81. https://doi.org/10.1007/s00170-021-07760-6.

[37] Trofimov A, Le Pavic J, Pautard S, Therriault D, Lévesque M. Experimentally validated modeling of the temperature distribution and the distortion during the fused filament fabrication process. Addit Manuf 2022;54:102693. https://doi.org/10.1016/j.addma.2022.102693.

[38] Zhang J, Van Hooreweder B, Ferraris E. T4F3: temperature for fused filament fabrication. Prog Addit Manuf 2022. https://doi.org/10.1007/s40964-022-00271-0.

[39] Owens JT, Das A, Bortner MJ. Accelerating heat transfer modeling in material extrusion additive manufacturing: from desktop to big area. Addit Manuf 2022;55: 102853. https://doi.org/10.1016/j.addma.2022.102853.

[40] Jayaraghul TK, Karthik K, Yaswanth A, Venkatesan M. Nozzle flow characteristics of P.E.E.K (Poly-ether ether ketone) material used in 3D-printing. Mater. Today Proc. 2021;44:2963–7. https://doi.org/10.1016/j.matpr.2021.02.004.

[41] Mishra AA, Momin A, Strano M, Rane K. Implementation of viscosity and density models for improved numerical analysis of melt flow dynamics in the nozzle during extrusion-based additive manufacturing. Prog. Addit. Manuf. 2022;7:41–54. https://doi.org/10.1007/s40964-021-00208-z.

[42] Serdeczny MP, Comminal R, Pedersen DB, Spangenberg J. Numerical simulations of the mesostructure formation in material extrusion additive manufacturing. Addit Manuf 2019;28:419–29. https://doi.org/10.1016/j.addma.2019.05.024.

[43] Tran TA, Lobov A, Kaasa TH, Bjelland M, Midling OT. CAD integrated automatic recognition of weld paths. Int J Adv Manuf Technol 2021;115:2145–59.

[44] Mokhtarian H. Product-process integrated meta-modeling using a graph-based approach. Tampere University; 2019.

[45] Mokhtarian H, Coatanéa E, Paris H, Mbow MM, Pourroy F, Marin PR, Vihinen J, Ellman A. A conceptual design and modeling framework for integrated additive manufacturing. J Mech Des 2018;140:1–13. https://doi.org/10.1115/1.4040163.

[46] Dua R, Rashad Z, Spears J, Dunn G, Maxwell M. Applications of 3d-printed peek via fused filament fabrication: a systematic review. Polymers (Basel) 2021;13:1–15. https://doi.org/10.3390/polym13224046.

[47] Verma S, Sharma N, Kango S, Sharma S. Developments of PEEK (Polyetheretherketone) as a biomedical material: a focused review. Eur Polym J 2021;147:110295. https://doi.org/10.1016/j.eurpolymj.2021.110295.

[48] Shekar RI, Kotresh T, Rao PD. Properties of high modulus PEEK yarns for aerospace applications. J Appl Polym Sci 2009;112:2658–67. https://doi.org/10.1002/app.29765.

[49] VICTREX. Victrex ® PEEK 450g. 2014.

[50] Dealy JM, Wang J. Melt rheology and its applications in the plastics industry, second. Springer Science 2013. https://doi.org/10.1007/978-94-007-6395-1.

[51] Kaseem M, Gun Y. Melt flow behavior and processability of polylactic acid / polystyrene (PLA / PS) polymer blends. J Polym Environ 2016. https://doi.org/10.1007/s10924-016-0873-5.

[52] Shibuya N, Porter RS. A kinetic study of PEEK sulfonation in concentrated sulfuric acid by ultraviolet- visible spectroscopy. Polymer (Guildf) 1994;35:3237–42.

[53] Bakrani Balani S, Chabert F, Nassiet V, Cantarel A. Influence of printing parameters on the stability of deposited beads in fused filament fabrication of poly(lactic) acid. Addit Manuf 2019;25:112–21. https://doi.org/10.1016/j.addma.2018.10.012.

[54] Escobedo J, Mansoori GA. Surface tension prediction for pure fluids. AIChE J. 1996;42:1425–33.

[55] Gibou F, Chen L, Nguyen D, Banerjee S. A level set based sharp interface method for the multiphase incompressible navier-stokes equations with phase change. J Comput Phys 2007;222:536–55. https://doi.org/10.1016/j.jcp.2006.07.035.

[56] Balani SBakrani, Cantarel A, Chabert F, Nassiet V. Influence of parameters controlling the extrusion step in Fused Filament Fabrication (FFF) process applied to polymers using numerical simulation. In: AIP Conf. Proc. 1960; 2018. https://doi.org/10.1063/1.5034995.

[57] COMSOL. Heat transfer module user's guide. 2018. https://doc.comsol.com/5.4/doc/com.comsol.help.heat/HeatTransferModuleUsersGuide.pdf.

[58] Bellehumeur CT, Kontopoulou M, Vlachopoulos J. The role of viscoelasticity in polymer sintering. 1998. https://doi.org/10.1007/s003970050114.

[59] Pokluda O, Bellehumeur CT, Vlachopoulos J. Modification of Frenkel's model for sintering. AIChE J. 1997;43:3253–6. https://doi.org/10.1002/aic.690431213.