Modelling Thermoplastic Filaments’ Sintering by Level Set Method

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): 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]: with α ∈ [0; 2.π[; R0, the final radius of the inverted ellipse (R0= √2.Ri, Ri being the initial filaments' radius [m]); parameter m ∈ [0; 1[, in Cartesian coordinates (Eq. 2), and parameters a (i.e., the size coefficient being equal to Lmax/2 with Lmax the maximal horizontal length of the sintered filaments; Fig.1-b), as well as B, the shape parameter [-] in polar coordinate system (Eq. 3); [rad], defined as α-π/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, tvs [s], from θ experimental values, as shown previously on a thermoplastic biopolymer [6,7]: , 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: � 2 � � � ℎ = ( 2 ). The assessment of tvs 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, tvs was found below 100s for temperature above 120°C. It is similar to values of tvs for ABS, an amorphous standard polymer, at its own printing temperature [6]. 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 twophase 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 (Øfilament ≈ 2mm). The amorphous extruded material, namely "Z20Gly", presents a glass transition temperature at Tg = 42°C, associated to a main mechanical relaxation at Tα ≈ 50°C, with material flowing above 80°C [6].
Hot melt sintering trials. Two adjacent filaments cut as cylinders (Lfilament = 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 -1 ; 8 bits coding; resolution about 200pixel.mm -1 ).
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 Lmax 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. 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: 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, ρ, the viscosity, μ, and the thermal diffusivity, D, at each point of the mesh are determined: (9) and (10) with ρ1 and ρ2, the densities; μ1 and μ2, the dynamic viscosities and D1 and D2, 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. = 0 for a mesh element filled with air (in blue, Fig.2-a) and = 1 for a mesh element filled with biopolymer, the interface being at = 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: with , the volume fraction [-]; γ, the re-initialization factor [m.s -1 ] (maximum speed rate of melt in TPF) and εIs, 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: with T, the temperature (Toven=120°C and Tinitial_Z20Gly=25°C); V, the flow velocity to link with NS equations; λ, the thermal conductivity [W.m -1 .K -1 ] and Cp, the specific thermal capacity [J.kg -1 .K -1 ].

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 2 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 ηZ20Gly (Table 1) from a factor (x1.0) to

Achievements and Trends in Material Forming
(x0.01), it increases the falling rate of the shape factor, m 2 , and decreases its final level, typical of a faster and more efficient sintering, respectively (Eq. 5, Fig.3-a). The slope in the linear decreasing domain of m 2 , as well as its final value in the plateau, is in a similar order of magnitude as the one obtained during sintering trials at 120°C with viscosity coefficient set at (x0.1). The characteristic viscous sintering time is then in the same range as the experimental one. With this setting, the evolution of the simulated images is found next to the experimental ones (Fig.3-c), but the simulated profiles show a contracting geometry, with 90% of the initial Lmax reached at the end of the simulated sequence. It is different from the experiment one, which presents an expansion of Lmax to about 115% of the initial value. The decrease of Lmax in the simulated conditions is due to the effect of the surface tension, the driving force of the sintering phenomenon, which brings together the two filaments as converging to a central position. In the experimental trial, the same phenomenon happens, but, probably due to viscoelastic swelling, leading to edge effects when monitored as a side view, and to gravity, with vertical melt flowing, a radial growth is globally observed on the acquired images [9].
The impact of the increase of the simulated environment temperature, from 100°C to 140°C, is found in a similar way as the one of the decrease of the viscosity, typical of its thermal dependency ( Table 1). The sensitivity study of the 2D-FEM model shows that at high temperature the sintering kinetics is enhanced, with an increase of the falling rate of m 2 (Fig.4-a). Results are also similar when comparing the evolution of simulated and experimental Lmax, with the effect of the surface tension reducing the radial size of the simulated sintered filaments: Lmax decreases to about 90% of its initial value after a tens seconds for simulated temperatures ranging from 100°C to 140°C (Fig.4-b). For the highest temperature, it has to be noticed that by the effect of the gravity and the vertical associated melt down-flowing increasing radial filaments' diameter, there is a slight initial increase of Lmax, up to 105%. Thus, the 2D modelling with constant amount of the two phases represented in one plan, does not allow reproducing neither the edge effect due to swelling nor impact of gravity on melt flowing. So, even if the whole sintering is simulated, such lacks may be at the root of the differences between experimental and simulated images [9]. 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-b- 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-).
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 Γ = 0.03N.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 Lmax, 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). a-b- Figure 6. Evolution of the parameter m 2 during experimental sintering at 120°C and simulations for various surface tension levels (a-). Evolution of Lmax, normalized as the ratio Lmax/Lmax_initial (b-).
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 tvs≈50s, obtained for similar tested parameters of viscosity and surface tension presented above. The impact of the decrease of viscosity, which enhances the fusionbonding phenomenon, is similar to the previous simulated one with non-refined meshing (Fig.3). a-b- 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 Lmax, normalized as the ratio: Lmax/Lmax_initial (b-).

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.