# Modelling Thermoplastic Filaments' Sintering by Level Set Method

Laurent Chaunier, Sofiane Belhabib, Sofiane Guessasma, Anne-Laure Reguerre, Eric Leroy

► To cite this version:

Laurent Chaunier, Sofiane Belhabib, Sofiane Guessasma, Anne-Laure Reguerre, Eric Leroy. Modelling Thermoplastic Filaments' Sintering by Level Set Method. Key Engineering Materials, 2022, 926, pp.365 - 372. 10.4028/p-710m7b. hal-03747252

HAL Id: hal-03747252

https://hal.science/hal-03747252v1

Submitted on 8 Aug 2022

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
# Modelling Thermoplastic Filaments' Sintering by Level Set Method

Laurent CHAUNIER¹, a*, Sofiane BELHABIB², b, Sofiane GUESSASMA¹, c,
Anne-Laure REGUERRE¹, d and Eric LEROY³, e

¹INRAE, BIA, UR 1268, F-44316, Nantes, France

²Université de Nantes, Oniris, CNRS, GEPEA, UMR 6144 F-44000 Nantes, France

³GEPEA, UMR CNRS 6144, F-44600 Saint Nazaire, France

a*laurent.chaunier@inrae.fr, bsofiane.belhabib@univ-nantes.fr, csofiane.guessasma@inrae.fr, danne-laure.reguerre@inrae.fr, eeric.leroy@univ-nantes.fr

**Keywords:** Biopolymer; FEM; Fusion-bonding; Morphological image analysis; Zein.

**Abstract.** The viscous sintering kinetics of thermoplastic polymers is generally studied by monitoring the evolution of the bonding neck between two particles (spherical, or cylindrical) and using a refined Frenkel-Eshelby's model. Recently, we showed that the entire contour of sintering filaments could be modelled by lemniscates as figure-eight shape curves to assess bonding abilities of a 3D-printable plasticized biopolymer. Using COMSOL Multiphysics® software, we set up a 2D finite element model of thermoplastic filaments' viscous sintering with flow front tracking by the level set method. This leads to contrasted images of the two phases, i.e. air and polymer, allowing the prediction of the shape of the interface corresponding to the filaments' contour. An image analysis procedure is applied to the simulated sequences and the ones acquired during sintering trials of extruded filaments based on zein, a corn protein plasticized by 20w% glycerol. This method is based on the assessment of the coordinates of sintered filaments' edge pixels and their fitting by lemniscates of Booth. We show that the 2D FEM approach combined with level set method allows simulating the hot melt viscous sintering of a 3D-printable thermoplastic biopolymer as a two-phase flow. Furthermore, the image analysis is successfully applied to simulated and experimental sequences, thanks to the monitoring of the filaments' contour, to assess their bonding kinetics and check its modelling.

## Introduction

Additive manufacturing techniques based on polymers' sintering as powders (e.g. Selective Laser Sintering, SLS) or filaments (e.g. Fused Filament Fabrication, FFF) lead to 3D printed objects obtained layer-by-layer thanks to molten material's bonding abilities. In the case of FFF, filaments' cross section may be considered as a representative basic element of the porous structure of 3D printed parts [1]. Bonding properties are at the root of the printability of a considered thermoplastic polymer and have an impact on the quality of the final object, as well as its end-use properties. The characterization of the viscous sintering is carried out thanks to the analysis of the length of the bonding neck, x [m], between two typical objects: cylindrical filaments or spherical geometries of diameter d [m]. The bonding angle θ [rad] is given by equation (1):

$$ \theta = \sin^{-1}\left(\frac{x}{d}\right) \qquad (1) $$

Based on such approach, a lemniscate model was proposed and solved numerically by Hopper [2] based on refined Frenkel-Eshelby's model [3,4] (Fig.1). But, the entire geometry of the system, i.e., the complete contour profile of two coalescing filaments, may be taken into account. Indeed, theoretical lemniscates of Booth can be fitted to the sintered filaments' entire contour. Such mathematical functions, leading to plane algebraic figure-eight shape curves, can be expressed in Cartesian coordinates and also be described in a polar coordinate system [5]:

$$ \begin{cases} x(\alpha) = R_0 \left[ (1 - m^2) \left( (1 + m^2)^{-\frac{1}{2}} \cos 2\alpha + m^2 \right)^{-1} \right] \cdot (1 + m) \sin\alpha \\ y(\alpha) = R_0 \left[ (1 - m^2) \left( (1 + m^2)^{-\frac{1}{2}} \cos 2\alpha + m^2 \right)^{-1} \right] \cdot (1 - m) \cos\alpha \end{cases} \quad \text{and} \quad r_p(\varphi) = a \cdot \sqrt{\cos^2 \varphi + B^2 \sin^2 \varphi} \quad (2) \text{ and } (3) $$
with $\alpha \in [0; 2\pi[$; $R_0$, the final radius of the inverted ellipse ($R_0=\sqrt{2}R_i$, $R_i$ being the initial filaments' radius [m]); parameter $m \in [0; 1[$, in Cartesian coordinates (Eq. 2), and parameters a (i.e., the size coefficient being equal to $L_{max}/2$ with $L_{max}$ the maximal horizontal length of the sintered filaments; Fig.1-b), as well as B, the shape parameter [-] in polar coordinate system (Eq. 3); $\varphi$ [rad], defined as $\alpha-\pi/2$ from Hopper's work, to match as an horizontal representation in polar coordinate system.

Such model allows the estimation of the characteristic time for viscous sintering, $t_{vs}$ [s], from $\theta$ experimental values, as shown previously on a thermoplastic biopolymer [6,7]:

$$t_{vs} = \left(\frac{R_i \cdot \mu}{\Gamma}\right) \quad \text{and} \quad \left(\frac{dm^2}{d\left(\frac{t}{t_{vs}}\right)}\right)_{\text{theoretical}} = \left(\frac{dm^2}{dt}\right)_{\text{experimental}} \cdot t_{vs} \quad (4) \text{ and } (5)$$

$t_{vs}$ is linked to material's thermo-rheological properties: $\mu$, the melt's viscosity [Pa.s] and $\Gamma$, its surface tension [N.m$^{-1}$]. This experimental slope, $(\frac{dm^2}{dt})_{experimental}$, is determined during sintering trials in an instrumented oven. The use of Hopper's abacus [2] makes it possible to obtain a theoretical master curve and its polynomial fitting: $(\frac{dm^2}{d\left(\frac{t}{t_{vs}}\right)}_{theoretical} = f(m^2)$. The assessment of $t_{vs}$ allowed defining the processing window of zein, a corn storage protein, plasticized by 20w% glycerol to obtain its 3D printing by FFF. In the case of zein, $t_{vs}$ was found below 100s for temperature above 120°C. It is similar to values of $t_{vs}$ for ABS, an amorphous standard polymer, at its own printing temperature [6].

a-

Extruded filaments

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

b-

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Scheme of the instrumented oven for images acquisition during isothermal sintering trials (a-). Dimensional characterization of filaments considering a side view and 4 lemniscates of Booth, as inverse ellipses mathematical functions, representing filaments' profile during theoretical hot melt sintering ($m = 0.95, 0.75, 0.50$, and $0.05$; b-).</figcaption></figure>

Modelling the hot melt sintering of a thermoplastic polymer has been investigated by Finite Element Method (FEM), especially in the case of circular powder particles [8]. Such 2D FEM model can be combined with the Level Set (LS) method to track the moving interface in the case of a two-phase flow, i.e. the molten biopolymer and air. The simulated sintering could be modelled as successive “8” shape profiles characterized by lemniscates of Booth fitted to their edges. Furthermore, the characteristic time for viscous sintering can be computed by the use of the Hopper's abacus to check the various settings of the 2D FEM model, as material's properties (e.g. viscosity and surface tension), processing parameters (e.g. temperature) and mesh size. Then, the aim of the present work is to investigate such modelling in the case of the isothermal sintering of biopolymer-based extruded filaments, thanks to the image analysis of their evolving geometry as successive lemniscates of Booth.
Experimental and Model Settings

**Biopolymer-based material.** A powder blend based on zein with 20w% glycerol, added as plasticizer, is extruded thanks to a twin-screw micro-compounder at 130°C (Haake Minilab, Thermo Scientific GmbH; Karlsruhe, Germany) through a cylindrical die (Ø<sub>filament</sub> ≈ 2mm). The amorphous extruded material, namely “Z20Gly”, presents a glass transition temperature at T<sub>g</sub> = 42°C, associated to a main mechanical relaxation at T<sub>α</sub> ≈ 50°C, with material flowing above 80°C [6].

**Hot melt sintering trials.** Two adjacent filaments cut as cylinders (L<sub>filament</sub> = 5mm), are placed in an instrumented oven with a transparent glass-made front. The isothermal viscous sintering is followed up thanks to a CMOS camera (1 frame.s<sup>-1</sup>; 8 bits coding; resolution about 200pixel.mm<sup>-1</sup>).

**Image analysis.** The same procedure is applied on the acquired images during sintering trials [7] and on the ones obtained by 2D FEM model combined with a Level Set method. A segmentation is carried out with ImageJ (free software; National Institutes of Health, USA) to obtain binary images with one white object, corresponding to the filaments (pixels=1), on a black background, corresponding to the backlighting (pixels=0). On each binary image of sintering sequences, a morphological image analysis is carried out to determine the pixels' contour of the two sintered filaments and to fit lemniscates. This procedure is implemented on Matlab® (The MathWorks Inc., Natick, MA, USA): each image is vertically centred on the ordinate at which the horizontal length of the white object is maximum (where L<sub>max</sub> is determined) and, horizontally, at the abscissa of the centroid, equivalent to a central point for the two coalescing filaments (Fig.1-b). In the case of a sole filament, the morphological characterization is carried out on binary images obtained thanks to the same instrumented oven.

**2D FEM model combined with Level Set method.** Numerical simulation of filaments' sintering is performed using two-phase flow simulation (TPF) with COMSOL Multiphysics® software (COMSOL Inc., Burlington, MA, US). The phase#1 is air and phase#2 is the thermoplastic biopolymer material (Fig.2). The software makes it possible to solve systems of differential equations with the coupling of physical phenomena by the Finite Element Method (FEM). The integrated solver allows the resolution of the equations and the visualization of the results by a graphic post-processing.

a-

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

b-

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Geometry and description of the domain for the numerical modelling of filaments' sintering as a two phase flow (Ø<sub>disk Z20Gly</sub> = 2 mm, each) (a-). Adaptive meshing of the modelled domain (b-).</figcaption></figure>

The level set (LS) and Navier-Stokes (NS) equations are used to model the flow of the molten biopolymer during hot melt sintering under the effect of the surface tension as the driving force. The monitoring of the interface between the two immiscible fluids is carried out between the biopolymer (phase#2), flowing at temperature above 80°C [6], and air (phase# 1). A reduced form can be given to the NS equation, because of inertial term being negligible at low Reynolds number:

$$ \rho \cdot \frac{\partial v}{\partial t} + \rho(V \cdot \nabla v) = \nabla \cdot [-pI + \mu(\nabla v + (\nabla v)^T)] + \rho \cdot g + F_{ST} + F \quad (6) $$

with t, the time [s]; ρ, the density [kg.m<sup>-3</sup>]; V, the flow velocity [m.s<sup>-1</sup>]; P, the pressure applied to the fluid [Pa]; μ, its dynamic viscosity [Pa.s]; g, the acceleration of the gravity field [9.80665 m.s<sup>-2</sup>, or N.kg<sup>-1</sup>], F<sub>ST</sub> the force resulting from the surface tension [N] and F the whole other forces applied to the flowing fluid [N]. The conservation of mass is solved by the continuity equation:
$$ \rho(V \cdot \nabla V) = 0 \qquad (7) $$

NS and continuity equations are solved during the simulation of the sintering of the two extruded filaments, applying the principles of conservation of momentum and mass, respectively. The density, $\rho$, the viscosity, $\mu$, and the thermal diffusivity, $D$, at each point of the mesh are determined:

$$ \rho = (1 - \varphi) \cdot \rho_1 + \varphi \cdot \rho_2, \quad \mu = (1 - \varphi) \cdot \mu_1 + \varphi \cdot \mu_2 \quad \text{and} \quad D = (1 - \varphi) \cdot D_1 + \varphi \cdot D_2 \quad (8), (9) \text{ and } (10) $$

with $\rho_1$ and $\rho_2$, the densities; $\mu_1$ and $\mu_2$, the dynamic viscosities and $D_1$ and $D_2$, the thermal diffusivities of phase #1 and #2, respectively.

The volume fraction is one of the outputs of the two-phase flow simulation to follow the simulated displacement (sintering) of the fluid volume and the air/melt interphase. $\varphi = 0$ for a mesh element filled with air (in blue, Fig.2-a) and $\varphi = 1$ for a mesh element filled with biopolymer, the interface being at $\varphi = 0.5$. The melt flow velocity is considered zero when the cylinders are in contact with the wall of the furnace on which they are placed. The monitoring of the interface between the two phases is carried out by the LS method by solving the following transport equation added to those of NS. It allows describing fluid flow in a continuous medium with conservation of the mass:

$$ \frac{\partial \varphi}{\partial t} + \nabla \cdot (V \varphi) = \gamma \nabla \cdot (\varepsilon_{ls} \nabla \varphi - \varphi (1 - \varphi) \frac{\nabla \varphi}{|\nabla \varphi|}) \qquad (11) $$

with $\varphi$, the volume fraction [-]; $\gamma$, the re-initialization factor [m.s⁻¹] (maximum speed rate of melt in TPF) and $\varepsilon_{ls}$, regulating the interfacial thickness (reduced for a clear separation of the two phases).

To take into account heat transfers during the simulation, allowing the calculation of the evolution of the temperature of the material in the model, a convection-diffusion equation is added to the model:

$$ \rho \cdot C_p \cdot \left( \frac{\partial T}{\partial t} + V \cdot \nabla T \right) = \nabla \cdot (\lambda \nabla T) \qquad (12) $$

with T, the temperature ($T_{oven}$=120°C and $T_{initial\_Z20Gly}$=25°C); V, the flow velocity to link with NS equations; $\lambda$, the thermal conductivity [W.m⁻¹.K⁻¹] and $C_p$, the specific thermal capacity [J.kg⁻¹.K⁻¹].

Concerning air (i.e., Phase#1), its thermal conductivity is assumed at $\lambda_{Air}$=0.024W.m⁻¹.K⁻¹, $\rho_{Air}$=1.2kg.m⁻³, $C_{p\_Air}$=1kJ.kg⁻¹.K⁻¹ and $\mu_{Air}$=1.10⁻⁵Pa.s. On the other hand, the density of the plasticized zein (i.e., Phase#2) is set at $\rho_{Z20Gly}$=1200 kg.m⁻³. $C_{p\_Z20Gly}$, previously determined by DSC, increases when exceeding the glass transition temperature from 1.5.10³ J.kg⁻¹.K⁻¹, at 20°C, to 2.2.10³ J.kg⁻¹.K⁻¹, at 130°C. The thermal conductivity is determined experimentally (Hotdisk® system) at $\lambda_{Z20Gly}$=0.26W.m⁻¹.K⁻¹, for a diffusivity at $D_{Z20Gly}$ = 1.34.10⁻⁷ m².s⁻¹. The measurement of melts' viscosity being tricky, due to proteins aggregation for long time at high temperature, it is deduced from characteristic viscous sintering time, $t_{vs}$, and according to Hopper's model [6] (Table 1).

<table><caption>**Table 1.** Characteristic viscous sintering time, $t_{vs}$ [s], and viscosity, $\eta$ [Pa.s], evaluated from sintering trials carried out at various temperatures for the plasticized zein according to Hopper's approach [6] (with $\eta=\Gamma \cdot t_{vs}/a$; with a, the radius of the filament [m] and $\Gamma$ assumed to be equal to 3.10⁻² N.m⁻¹).</caption><thead><tr><th>Temperature [°C]</th><th>t<sub>vs</sub> [s]</th><th>η<sub>Z20Gly</sub> [Pa.s]</th></tr></thead><tbody><tr><td>80</td><td>3492</td><td>1,1.10<sup>5</sup></td></tr><tr><td>100</td><td>380</td><td>1,1.10<sup>4</sup></td></tr><tr><td>120</td><td>82</td><td>2,5.10<sup>3</sup></td></tr><tr><td>125</td><td>30</td><td>8,9.10<sup>2</sup></td></tr><tr><td>130</td><td>25</td><td>7,5.10<sup>2</sup></td></tr><tr><td>140</td><td>14</td><td>4,3.10<sup>2</sup></td></tr></tbody></table>

## Results and Discussion

The numerical simulation leads to images with neat interface between the thermoplastic material and the surrounding air, allowing the computation of m² and Lmax parameters of the fitted lemniscates during simulated sintering sequences, as previously done on experimental trials (Fig.3). In the case of a sensitivity study of the 2D-FEM model, the viscosity is shown to greatly influence sintering kinetics. Indeed, by decreasing the coefficient applied to $\eta_{Z20Gly}$ (Table 1) from a factor (x1.0) to
a-

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

b-

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of the parameter $m^2$ during filaments’ experimental sintering at 120°C and simulations for various temperatures (a-). Evolution of $L_{max}$ normalized as the ratio: $L_{max}/L_{max\_initial}$ (b-).</figcaption></figure>

To go further in the modelling of filaments sintering, the evolution of the shape of a sole filament has been reproduced by a numerical simulation. The vertical flowing due to gravity is well reproduced in this case, where coalescence does not occur (Fig.5-a). During experimental trials, the average perimeter of a sole filament is increased up to 115% of its initial value, even if standard deviation are non-negligible (Fig.5-b). In the case of the 2D-FEM model, the perimeter remains stable. Such difference is due to edge effects, at the root of the increase of filaments’ radial size in the experimental trials, which is not possible in the simulation as 2D side-view without adding material in the 2D plan. As previously found during sintering sequences in the case of two adjacent filaments, the viscoelastic swelling due to 3D displacement of the molten material, from the middle of the filament to its edges, leads to the swelling of a sole filament placed at high temperature when experimentally followed up as a side-view imaging.

a-

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

b-

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Results of the simulation for a sole filament at 120°C as successive binary images (a-). Evolution of the perimeter of a sole filament placed at 120°C in the instrumented oven and simulations for various viscosity coefficients in the case of a sensitivity study (b-).</figcaption></figure>

The sensitivity study of the 2D FEM model combined with LS method has been extended to the variation of the surface tension, from a factor (x0.5) to (x2), based on a basic value set at $\Gamma = 0.03\text{N}\cdot\text{m}^{-1}$ (Fig.6). By limiting the value of such a driving force of the sintering phenomenon, in the case of factor (x0.5), the driving force for sintering in lowered, decreasing the falling rate of $m^2$ and the vertical melt flow due to gravity (Fig.6-a and 6-b, respectively). It is combined due to the effect of gravity that is favoured as such low surface tension level, leading to increase the radial (horizontal) size of the sintered filaments with an increase up to 106% of the initial filaments’ diameter (Fig.6-b). Then the coalescence leads to a decrease of $L_{max}$, because of the displacement of the filaments one to the other, as a contracting size of the modelled molten material. Such size reduction, down to 75% of the initial size after 4 minutes sintering, is enhanced when increasing the surface tension up to the maximum value of the applied factor (x2).
<figure><img src="image_10.png" /><figcaption>b-</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Evolution of the parameter $m^2$ during experimental sintering at 120°C and simulations for various surface tension levels (a-). Evolution of $L_{max}$, normalized as the ratio $L_{max}/L_{max\_initial}$ (b-).</figcaption></figure>

Finally, the numerical simulation has been also carried out with a refined meshing of the modelled domain, after dividing the mesh size by a factor (/2) (Fig.7). The kinetics is less jagged concerning the evolution of $m^2$ lemniscates’ model parameter, but similar results are obtained in term of its decreasing slope, with the same $t_{vs} \approx 50$s, obtained for similar tested parameters of viscosity and surface tension presented above. The impact of the decrease of viscosity, which enhances the fusion-bonding phenomenon, is similar to the previous simulated one with non-refined meshing (Fig.3).

<figure><img src="image_12.png" /><figcaption>b-</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Evolution of the parameter $m^2$ during experimental sintering at 120°C and simulations with refined meshing for various viscosity level (a-). Evolution of $L_{max}$, normalized as the ratio: $L_{max}/L_{max\_initial}$ (b-).</figcaption></figure>

## Conclusion

The image analysis procedure with lemniscates fitting is efficient to characterize both filaments’ shape and size factors during hot melt sintering sequences from experimental trials and numerically simulated by 2D FEM model combined with the Level Set method. A sensitivity study of the model shows that the decrease of melts’ viscosity, as well as the increase of their surface tension, enhance sintering kinetics, but also that 2D approach presents a lack in considering melts viscoelastic swelling and its associated edge effects. This phenomenon affects the characterization of the fusion-bonding abilities, as well as its modelling. Thus, the 3D monitoring of the whole filaments’ envelope has to be carried out, per example by optical profilometry or X-ray computed tomography coupled to heating systems, and then implemented in hot melt sintering models for better matching thermoplastic materials’ behaviour.
References

[1] B.N. Turner, S.A. Gold, A review of melt extrusion additive manufacturing processes: II. Materials, accuracy, and surface roughness. Rapid Prototyping Journal 21 (2015) 250-261.

[2] R.W., Hopper, Coalescence of Two Equal Cylinders: Results for Creeping Viscous Plane Flow Driven by Capillarity. Com. Am. Cer. Soc. 67 (1984) C262-C264.

[3] C. Bellehumeur, L. Li et al., Modeling of bond formation between polymer filaments in the fused deposition modeling process, J. Man. Proc. 6 (2004) 170-178.

[4] N.D. Polychronopoulos, J. Vlachopoulos, The role of heating and cooling in viscous sintering of pairs of spheres and pairs of cylinders, Rap. Prot. J. 26(4) (2020) 719-726.

[5] K. Sowinski, P. Jasion, Strength and stability of shells based on Booth lemniscate loaded with external pressure. Thin-Wall. Struc. (2019) 144.

[6] L. Chaunier, G. Della Valle et al., Viscous sintering kinetics of biopolymer filaments extruded for 3D printing. Polym. Test., (2019) 77 105873.

[7] L. Chaunier, A.-L. Réguerre, E. Leroy, Characterization of thermoplastic filaments sintering kinetics by image analysis and adjustment of Booth lemniscates, ESAFORM Proceedings (2021).

[8] C. Balemans, M.A. Hulsen, P.D. Anderson, Sintering of Two Viscoelastic Particles: A Computational Approach, Appli. Sci. (2017) 7(5) 20.

[9] S. Bordere, D., Gendron et al., Monte Carlo prediction of non-Newtonian viscous sintering: Experimental validation for the two-glass-cylinder system. J. Am. Cer. Soc. (2005) 88(8), 2071-2078.