Next Article in Journal
Valorization of Banana Peel Using Carbonization: Potential Use in the Sustainable Manufacturing of Flexible Supercapacitors
Previous Article in Journal
Single-Shot Multi-Frequency 3D Shape Measurement for Discontinuous Surface Object Based on Deep Learning
Previous Article in Special Issue
High Precision 3D Printing for Micro to Nano Scale Biomedical and Electronic Devices
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Experimental and Numerical Investigation of the Die Swell in 3D Printing Processes

Dipartimento di Ingegneria Chimica, dei Materiali e Della Produzione Industriale, Università Degli Studi di Napoli Federico II, P.le Tecchio 80, 80125 Napoli, Italy
*
Author to whom correspondence should be addressed.
Micromachines 2023, 14(2), 329; https://doi.org/10.3390/mi14020329
Submission received: 11 January 2023 / Revised: 20 January 2023 / Accepted: 21 January 2023 / Published: 27 January 2023
(This article belongs to the Special Issue 3D Printed Micro-/Nano Devices)

Abstract

:
Fused deposition modelling is one of the most widely used additive manufacturing techniques and the diffusion of 3D printers has increased in popularity even further in recent times. Since high precision is required in 3D printing, a good control over the extrusion process is necessary. In this regard, a crucial phenomenon to be accounted for is the die or extrudate swell, i.e., the enlargement of the cross-section of the strand when coming out of the printer nozzle. While this phenomenon has been studied in large scale extruders, it has not yet been investigated in depth for 3D printing processes. In this work, the die swell phenomenon observed in a printed PLA filament is studied by experiments and fluid dynamic simulations. A novel, easy-to-use, accurate and fast procedure for measuring the value of the die swell ratio during the printing process is developed, accounting for typical errors related to a non-constant strand diameter and possible oscillations of the filament with respect to the extrusion direction. As the printing velocity is increased, a linearly increasing swelling ratio is observed at low printing speeds. The trend flattens at moderate speed values. A sudden increase is found at high printing velocities. The swelling ratio measured with the proposed technique is compared with the results of multi-mode viscoelastic simulations at different temperatures. A fair agreement between the experimental measurements and the numerical predictions is found for printing velocities that are typically employed in commercial 3D printers, supporting the reliability of the developed procedure.

1. Introduction

Additive manufacturing (AM) techniques have been employed both at industrial level from medium to high scales, and at lower scales in households, offices, and research centers. As a result of their widespread use in a variety of applications, the market for AM has been rapidly growing in recent years [1,2]. In this context, one of the most widely used techniques is fused deposition modelling (FDM), also known as material extrusion or fused filament fabrication. This technique builds 3D objects starting from a solid polymeric filament that is fed to a heated liquefier, acting both as the feedstock for the extrusion process and as a piston that provides the pressure allowing the material to flow. While the filament is being melted by the heat provided in the liquefier, it is forced through a nozzle and extruded. Three-dimensional printers are machines commonly used in material extrusion in which a motor moves the print nozzle in the three spatial directions as the strand of material is deposited on a flat surface. Three-dimensional objects of many complex shapes are then built by deposition of further layers that solidify on top of each other. For this reason, the filaments are usually made of amorphous or slow crystallising polymers, since fast crystallization would not allow the freshly deposited strand to weld with the previously printed layer [3]. The most commonly used polymer in FDM is polylactic acid (PLA), on which many research works have focused on [4,5].
The reduction in costs and control improvements for extruders and 3D printers have strongly contributed to the increase in popularity of FDM techniques in recent years. The future growth of FDM is, however, still limited by the understanding of the extrusion and deposition processes and the improvement in the manufacturing systems [1,6]. The complexity of the large number of physical phenomena involved requires a strong research effort to improve the performance and optimize the process. In this framework, both experimental and simulation works are useful.
One crucial aspect of controlling the extrusion processes is the understanding of the flow behavior both inside and outside the liquefier, linked to all the properties of the viscoelastic fluids employed in 3D printing that have a nonlinear dependence on shear rate and temperature. Once the filament is melted, the fluid encounters a contraction in the conical part of the nozzle and undergoes both shear and elongational stresses [3]. When the polymer melt is released from the nozzle, it is no longer constrained by the wall and the presence of a free surface allows it to release its stored elastic energy by enlarging its transversal size, giving rise to a phenomenon known as die swell or extrudate swell. In a Newtonian fluid, the die swell is explained by the rearrangement of the velocity profile, while in viscoelastic fluids, the extrudate swell is mainly attributed to elastic recovery [7,8,9].
A current challenge in FDM processes is the design of the deposition strategy, taking into consideration the polymer extrudate swell that affects the size of deposited strands and the accuracy of the built structure [10,11,12]. When the swell is not considered, the quality of the 3D printed objects is poor [13]. Despite its relevance, the extrudate swell has rarely been investigated in 3D printers. Currently, the design of the deposition strategy relies only on empirical methods and the effect of the processing conditions (e.g., shear rate and temperature) is neglected in standard FDM software. In addition, despite the many investigations that have been carried out on flow defects, the measurements of die swell are often not consistent among different experimental works [14,15,16]. This is due to the absence of a standard procedure for measuring the equilibrium value of the diameter because of a transient zone between the nozzle and the equilibrium of the strand and the existence of oscillations in the strand and in the temperature of the nozzle. It is worth noting that, when increasing the flow rate, three kinds of flow instabilities may occur. First, small defects on the extrudate surface, known as sharkskin, are encountered. Then, at higher flow rates, pressure oscillations cause stronger surface defects. Finally, if the flow rate is further increased, melt fractures occur [17]. Of course, the onset of melt fractures sets an upper limit for swelling measurements.
Several simulations regarding the die swell in extrusion processes have been performed, investigating different aspects of the phenomenon. The first numerical works date back to over 30 years ago and were limited to 2D geometries and simple constitutive equations due to the reduced computational power available at that time [18,19,20]. In recent years, extrudate swell simulations have become increasingly more complex, accounting for accurate viscoelastic and viscoplastic models [21,22], 3D channel geometries [23,24,25], temperature variation, and thixotropy [26]. It is broadly acknowledged that single mode viscoelastic models do not describe accurately enough the non-linear properties of polymeric melts [27]; thus, a quantitative prediction of the extrudate swell phenomenon requires the use of multi-mode viscoelastic constitutive Equations [28,29,30]. Only a few works, however, have compared experimental evidence with simulation results [28,30]. Furthermore, to the authors best knowledge, the analyses of the die swell in commercial 3D printers at small characteristic scales are much more limited.
In this work, we carried out a detailed analysis of the die swell from PLA extrusion footage in a commercial 3D printing device through both experiments and numerical simulations. We developed a novel procedure to accurately measure the diameter of the extruded strand in time at different printing velocities and temperatures, accounting for all the possible issues that might affect the accuracy of the measurements. We have compared the experimental outcomes with fluid dynamic simulation results by employing a multimode viscoelastic constitutive equation with parameters obtained by fitting the PLA rheology.

2. Materials and Methods

2.1. Experiments

2.1.1. Materials

Polylactic acid (PLA), grade 710M, was supplied by BewiSynbra (Sweden) in the form of pellets. This grade is characterized by a content of the isomer D-PLA equal to ca. 5.5 wt% and a polydispersity index (PI) equal to ca. 3. The PLA pellets were dried overnight at 60 °C in an oven under vacuum conditions before any manipulation. The drying protocol was based on work in the literature data that reported an efficient removal of humidity after 12 h at 60 °C [31]. A PLA filament with a diameter of 1.75 mm was produced using a Composer350 extruder (from 3DEVO company, Utrecht, The Netherlands). The diameter of the filament was measured and controlled by the Composer350 extruder setup.

2.1.2. Experimental Apparatus

A Prusa Research I3 Mk3s 3D printer was used to extrude the filament, shown in Figure 1. The experimental setup consists of a commercial 3D printer in which three zones can be identified: (a) the feeding zone, where the PLA filament is fed at constant speed by setting the speed of two counter-rotating wheels through a stepper motor (Nema 17) (Stepperonline, New York, NY, USA); (b) the melting zone, where the filament is fused by heat conduction with the metal walls kept at a constant temperature with the action of a 12–140 W cylindrical resistor (RepRap, Wateringen, The Netherlands) and a NPT sensor (RepRap, Wateringen, The Netherlands) with an accuracy of one degree Celsius; and (c) the exit zone, where the extruded strand swells and is recorded by a CCD videocamera (DMX TIS 1/1.8″ CMOS 3072 × 2048 Monocro from DB Electronics, Milan, Italy) equipped with an optical lens (2.0×, 2/3″ SilverTL Telecentric Lens from Edmund, York, UK) aligned at the nozzle exit to record the fluid flow. The use of the 2.0× optical lens allows a resolution of 3 microns/pixel for the calculation of the volumetric flow rate. The background is illuminated by an LED diffused light.
The experiments were performed using a nozzle shaped as shown in Figure 2. The inlet diameter of the nozzle was D in = 2 mm, the outlet diameter was D die = 0.6 mm, the convergence angle, θ , was 60 ° and the land length, L die , was 1.2 mm. The printer was controlled by using a software called Pronterface (RepRap, USA) that allowed setting of the extrusion speed of the filament (in mm/min), the temperature of the metal nozzle (in °C) and the length to extrude (in mm). The images of the extruded strands at the nozzle exit were acquired by means of the Imaging Source software (The Imaging Source Europe GmbH, Bremen, Germany) at 17.3 frames per second which were saved as movies in AVI format.

2.1.3. Experimental Protocol

The experimental procedure consists of: (a) filament loading, (b) temperature rise and stabilization for at least 10 min after the desired value is reached and (c) extrusion at a constant extrusion speed for 100 mm and, in parallel, optical observation of the die swell through camera recording. We have chosen a design of experiments (DoE) in which the temperature, T, was set at three different values above the melting point of PLA, i.e., 160 °C, 180 °C and 200 °C, and the extrusion speed, u in , was varied by choosing 16 values from 5 to 500 mm/min. Both the temperature and extrusion speeds correspond to typical values in 3D printing processes. In total, the presented DoE includes a set of 48 different experiments and a subset of five experiments, randomly chosen, that were repeated to check the consistency of the results.
A consistent measurement of the die swell needs to consider the following aspects: (a) the effect of the gravitational forces on the shape of the strand, (b) the transient zone of the strand diameter from the nozzle exit to the steady zone, (c) the vibrations and oscillations of the exiting strand and (d) the oscillation of the flow rate and temperature due to the 3D printing process. In our experiments, the effect of the gravity was neglected because the viscous forces were higher than the gravitational ones, i.e., the Galilei number, G a = g L 3 / ν (where g, L and ν are the gravitational acceleration, the characteristic size of the strand and the kinematic viscosity, respectively), is much smaller than one.
We have developed a custom-made Matlab script to evaluate whether the strand diameter was constant (after the transient zone) and to provide a consistent measure of this diameter considering both the oscillations of the strand and any possible oscillation in time of the flow rate and temperature. The procedure performs the following operations to analyze the images: (i) It loads the recorded movie of the die swell and extracts all the images, as shown in Figure 3a. (ii) Each image is binarized using the function i m b i n a r i z e with a constant threshold value equal to 0.4; then, the operator needs to select a rectangle to crop the image in the zone where the diameter does not change qualitatively, as shown in Figure 3b. (iii) The profile of the strand is recognized and stored as two vectors, by means of the function b w b o u n d a r i e s , represented by the red lines in Figure 3c. (iv) Two diameters are measured on the cropped image: the horizontal diameter ( D h ) equal to the horizontal distance between two points with the same ordinate on the recognized red profile (shown in Figure 3c with a dashed green line) and the perpendicular diameter ( D p ), equal to the minimum distance between a point and the other side of the red profile (i.e., the transversal line to the strand axis, shown in Figure 3c with a dashed yellow line). D p is calculated as the product of the D h and the cosine of the angle formed between D h and D p . Hence, D p is equal to D h when the strand is perfectly vertical and it is smaller when the strand has a slope with respect to the extrusion direction. By calculating D p , it is possible to avoid any error in the diameter measurement due to strand oscillations.
For each image, the script calculated N values of D p and D h , where N is the vertical dimension, in pixels, of the cropped image. In our experiments, N ranged from 500 to 1000 depending on the manual cropping. These values were averaged and the standard deviation was calculated. A high standard deviation is either an indication that the surface is very rough (e.g., melt fracture) or that the strand is still in the transient region where the diameter is growing.

2.1.4. Fluid Rheology

The viscosity of PLA was measured with a Physica MCR 301 rheometer from Anton Paar, Germany. The measurements were performed under a nitrogen atmosphere at different temperatures (160 °C, 180 °C and 200 °C) using a disposable plate–plate system with a radius of 12.5 mm (PP25) under oscillatory conditions. The gap between the plates was set to 1 mm. The sample was changed for each measurement and at each temperature, and the complex viscosity was determined through a frequency sweep test after 10 min of waiting time to allow the sample temperature to stabilize. The strain amplitude was 10 % , and the angular frequency varied from 0.1 to 600 rad/s. The strain amplitude was determined by means of a dynamic strain amplitude sweep in order to be in the linear viscoelastic regime. The linear viscoelastic moduli, G and G , and the complex viscosity, η , are shown as symbols in Figure 4 for the three temperatures.

2.2. Numerical Simulations

2.2.1. Mathematical Model

Since the nozzle geometry has an axis of symmetry, the problem was simulated in a 2D axisymmetric domain, as shown in Figure 5. The thicker lines represent the boundaries of the nozzle, whereas the thinner ones are the boundaries of the external domain.
The dimensions of the domain were chosen to match the nozzle of the 3D printer and are reported in Table 1. The values of L in , L air and R air were chosen to be large enough to neglect any effect of the inlet and outlet boundary conditions on the simulation results. The origin of the axes was placed at the exit of the nozzle on the symmetry axis.
The simulations were carried out by solving the mass and momentum balance equations governing the fluid dynamics of the melted polymer through the nozzle and in the surrounding environment (air). Hence, a multiphase system (PLA + air) was considered. We assume incompressible fluid and laminar flow conditions, justified by the very low values of the Reynolds numbers (see later). Gravity effects were negligible as discussed above. Even though temperature variations were present in the problem, for the sake of simplicity, the fluid was assumed to be at the temperature imposed by the heated walls of the liquefier, thus the problem was considered isothermal. This is a reasonable assumption at relatively low flow rates, as the temperature in the liquefier has time to become homogeneous.
Under these assumptions, the mass and momentum balance equations for the PLA read as:
· u = 0
ρ u t + u · u = p + · σ
where u , p, ρ and σ are the fluid velocity, pressure, density and the viscoelastic stress tensor, respectively. The latter is expressed through sum of M modes:
σ = n = 1 M σ i
where the stress tensor of each mode is given by the Giesekus constitutive Equation [32,33]:
λ i σ i + σ i + α i λ i η i σ i · σ i = 2 η i D
In this equation, λ i , α i and η i are the relaxation time, the mobility and the polymer viscosity of the i-th mode, respectively, D is the rate-of-deformation tensor:
D = 1 2 ( u + ( u ) T )
and σ i is the upper-convected time derivative that assures frame invariance:
σ i = σ i t + u · σ i ( u ) T · σ i σ i · u
The use of a multi-mode Giesekus model to describe viscoelastic fluids is well documented in the literature for various flow problems [34,35,36].
The number of modes, M, and the constitutive parameters were obtained by fitting the model predictions with the rheological data presented above. Specifically, a regression of the linear viscoelastic moduli gives M, the relaxation times and the polymer viscosities. The mobilities were then estimated by fitting the complex viscosity data. We found that four modes were sufficient to accurately describe the PLA rheology. The model predictions are shown as solid curves in Figure 4. The values of the constitutive parameters at the three temperatures are listed in Table 2.
The equations governing the fluid dynamics of the air are the same as Equations (1) and (2) with σ = 2 η air D and ρ replaced by the air density, ρ air . The surface tension force, γ , acts at the liquid–air interface. The other physical parameters used in the simulations are listed in Table 3.

2.2.2. Boundary and Initial Conditions

With reference to Figure 5, a fixed velocity profile is set in inlet, Γ in , with a constant velocity, u in , in the axial direction z, as imposed from the printer control software. On the outlet boundary, Γ out , and on the boundaries of the air domain, Γ air , a constant pressure equal to the ambient one is applied. Due to the fluid incompressibility, we can set such a pressure level equal to 0. At the nozzle walls, Γ wall , no-slip conditions are imposed, whereas on the symmetry axis, Γ sym , a symmetry condition is set. Finally, a boundary condition for the extra stress tensor, σ , is required in inflow. We set σ = 0 . The choice of L in assures that the velocity and stress fields fully develop in the inlet channel before entering the contraction.
The problem was initialized with velocity, pressure and stress identically null in the whole domain. In the initial conditions, the PLA occupied the whole nozzle domain up to 0.05 mm from the the end of the die, whereas air fills the rest of the domain, as shown in Figure 6, in which PLA is represented in red and air in blue.

2.2.3. Dimensionless Numbers

From the model parameters in Table 2, we can compute the relevant dimensionless numbers for the problem; namely, the Reynolds, Weissenberg and capillary numbers:
R e = ρ u in D in η 0 W i = λ u in D in C a = η 0 u in γ
The characteristic length was chosen as the inlet diameter D in = 2 mm, while the characteristic velocity was the average inlet velocity, u in . The maximum relaxation time was selected to define the Weissenberg number and the zero-shear viscosity, η 0 , was used in the definition of the Reynolds number. Table 4 reports the values of R e , W i and C a for the maximum average velocity considered in this work and for the three investigated temperatures. It can be readily observed that the Reynolds number was very low, as typically occurs for these systems. A Weissenberg number of up to about eight was obtained, denoting the relevance of elastic effects. The capillary number was quite high, meaning that the viscous forces are much stronger than the surface tension ones.

2.2.4. Simulation Software

All the simulations have been performed using RheoTool (version 5.0) [37,38], which is a toolbox of the open-source software OpenFOAM, specifically designed to deal with viscoelastic fluids. The solver used is rheoInterFoam, a transient solver for two-phase laminar, isothermal flows that combines the interFoam solver of OpenFOAM with the rheoFoam solver of RheoTool to simulate multiphase viscoelastic flows. This solver is based on the volume-of-fluid (VOF) approach to capture the interface between the two phases. In this method, the mass and momentum balance equations of the two fluids are combined in a single set of equations. A new variable, termed ‘indicator function’ or ‘color function’, denoted by α , is added to identify the regions of the domain occupied by one of the two phases. Specifically, the PLA and air domains are identified by α = 1 and α = 0 , respectively. The liquid–air interface corresponds to α = 0.5 . This variable is transported by the fluid velocity so a convection equation for the indicator function is added to the mathematical model. The physical properties in the combined Equations (e.g., density and viscosity) are weighted between the properties of the single phases through the indicator function. The indicator function requires initial and boundary conditions. The initial condition is set as shown in Figure 6, i.e., we set α = 1 in all the mesh cells occupied by PLA (red) and α = 0 elsewhere (blue). The boundary conditions are α = 1 on Γ in , since only the viscoelastic fluid enters the domain, a symmetry condition on Γ sym and a zero-gradient condition on the other boundaries.
The software solves the transient problem described by the combined mass and momentum balance equations, the viscoelastic constitutive equations, and the convective equation for the indicator function. We set a sufficiently large final time to let the velocity and stress fields as well as the extrudate profile attain a steady-state. The die swell ratio is then computed as the ratio between the diameter of the extrudate in the zone where it has an horizontal profile and the nozzle diameter D die .

2.2.5. Mesh and Mesh Convergence Study

The domain was discretized in quadrilateral elements. A mesh convergence study was performed on a test case at  T = 180  °C and  u in = 100  mm/min using four different mesh resolutions that have the same element distribution but a different level of refinement. The details of the meshes used are reported in Table 5.
In Figure 7, the PLA–air interface outside the nozzle (obtained as the contour corresponding to α = 0.5 ) is plotted against the z-coordinate for the four different meshes. Apart from the mesh M1 which is extremely coarse, the extrudate profiles corresponding to the finer meshes are fairly superimposed. By increasing the inlet velocity, we found that the mesh M3 always produces results equivalent with mesh M4. The former was then used in the present work. The complete mesh is shown in Figure 8a together with a zoomed-in image of the die and of the exit zone in Figure 8b.

3. Results and Discussion

3.1. Experimental Results

Figure 9 shows the horizontal and perpendicular diameters of the strand, D h and D p , as function of the vertical length of the cropped image as obtained by the custom-made script. The origin of the abscissas in Figure 9 is set at the top horizontal boundary of the cropped image (blue box). Hence, the values on the x-axis of the graphs correspond to the vertical distance from the top boundary of the cropped image. Both diameters are normalized by the die diameter D die = 2 R die . Panels (a) and (b) refer to two images taken at different instants of the same experiment, i.e., 160 °C and 100 mm/min. In Figure 9a, the two sets of data are very close over the entire length of the cropped image. Indeed, the strand shown in the inset is almost vertical and its borders (denoted by red lines) are almost parallel to the yellow dashed axis that represents the extrusion direction. Moreover, the data in the graph are horizontal, meaning that the strand diameter is constant in the cropped region. On the contrary, the blue and red dots in Figure 9b are close at small lengths but deviate with increasing the length. In this case, the strand shown in the inset is tilted with respect to the extrusion direction. The increase of D h is due to the slope of the strand and not to a real increase in the diameter. The perpendicular diameter is, in fact, the most accurate measurement of the die swell because it takes under consideration the slope of the strand and its oscillations during measurement.
Figure 10 reports the results of the die swell as a function of the time for the experiment performed at 160 °C and 100 mm/min. Each point in Figure 10 is the average of the horizontal and perpendicular diameters over the length of the cropped image (as shown in Figure 9) and the shaded areas are the standard deviations. Notice that D h (blue dots) has a bigger error and higher fluctuation frequency with respect to D p (red dots) due to the aforementioned vertical oscillations of the strands during the experiments. However, the measurements show some fluctuations in D p as well. These can be attributed to the temperature fluctuation (with a period of around 10 s) in the process, related to the PID control of the temperature, as also reported in the literature [4]. The consistency of the measurements was checked by repeating the experiments three times at five different velocities for the lowest and the highest temperature. The results fall within the experimental error, shown in Figure 10.
In the rest of the work, the experimental results will be reported only in terms of D p as, as previously mentioned, it is a more accurate estimation of the real strand diameter.

3.2. Simulations Results

Simulations were performed with different values of inlet velocity ranging from 20 to 400 mm/min. The effect of the temperature and inlet velocity on the fluid rheological properties was evaluated by comparing, in dimensionless form, the relevant components of the velocity and stress tensor in the converging zone of the nozzle, in the die and in the exit zone of the fluid for three sets of T and u in .
As discussed in Section 2.2.3, the velocity components were made dimensionless by the inlet velocity and the stress components by u in η 0 / D in . In Figure 11, the dimensionless axial and radial components, u z and u r , of the velocity are shown. The differences inside the inlet channel and the first part of the die are minimal in all three simulations and for both axial and radial components, as expected for a fully developed flow. In the exit zone of the die and in the first part of the extruded strand, there are quantitative differences in the velocity. When the inlet velocity increases or the temperature decreases, the swelling of the viscoelastic fluid increases and, as a consequence, the radial component of the velocity increases in proximity of the channel exit. The axial component, instead, decreases in magnitude with the increase in the swelling ratio, and hence of the transversal area, due to the mass conservation. In all the simulations, the velocity profile becomes a plug-flow far from the die exit, as also reported in the literature [39].
In Figure 12, the relevant components of the stress tensor, σ zz , σ rr and σ rz . are plotted in dimensionless form for the same three simulations. For all three components, as the temperature decreases or the velocity increases, the stress increases both in the nozzle contraction and at the exit of the die. A large stress region is, in particular, observed at the contact point between the polymer and air at the end of the die due to the geometric singularity.

3.3. Die Swell Ratio: Experiments vs. Simulations

Figure 13 shows the die swell ratio, D / D die , at the three temperatures as a function of the printing speed. The circles are the experimental measurements and the stars denote the simulation results. For all the investigated temperatures, the die swell ratio increases monotonically as a function of the printing speed, while it decreases by increasing the temperature from 160 °C to 200 °C due to the reduction in the normal stresses. By looking at the experimental measurements, the data show an inflection point at a velocity of about 200 mm/min. In particular, three regions can be identified: (i) Region I, from 0 to about 100 mm/min (with a range that slightly depends on the temperature), is characterized by an almost linear increase in the die swell ratio with a slope that decreases with increasing temperature. (ii) Region II, from 100 to about 350 mm/min, is characterized by a very weak increase in the die swell ratio with the printing velocity. (iii) Region III, from 350 to 500 mm/min, is characterized by a sudden increase in the die swell ratio. Notice that, in this latter region, the error bars are much larger than those obtained at lower velocities.
In the first two regions, a remarkable agreement between the experiments and the simulations is found. Specifically, the simulations were able to capture both the nearly linear increase in the swelling ratio observed in Region I as well as the flatter trend noticed in Region II. The agreement is also quantitative as most of the simulation points fall within the experimental error bars for all the investigated temperatures. On the other hand, the abrupt increasing trend measured in the experiments at higher velocities (Region III) is not captured by the simulations, as they predict a nearly constant die swell ratio, similar to that in Region II.
A possible justification of the observed results is that at higher printing velocities, the residence time of the filament in the liquefier is not sufficient to reach a homogeneous constant temperature in the entire cross-section of the strand (i.e., the residence time is lower than the heating time). Thus, the sudden increase in the die swell slope in Region III observed in the experiments is likely due to the strong temperature inhomogeneity within the filament, as also reported in the literature [40,41,42], making, in fact, the heat transfer a limiting factor in polymer extrusion, especially at high flow rates. Such a temperature inhomogeneity was not considered in the simulations. Furthermore, even if in the investigated range of printing velocity no melt fracture phenomena were observed, the minor flow instabilities that occur are not easily captured with numerical simulations. It is hence evident how the simulations are able to describe the printing process at low and moderate extrusion velocities; whereas, at higher flow rates, the isothermal assumption is no longer suitable to accurately describe the process. Finally, it has to be mentioned that the die swell phenomenon depends on the shear as well as extensional stresses to which the fluid is subjected inside the nozzle [7]. The calibration of the constitutive parameters was, however, only based on shear viscosity data which is extensional data that is not available and much more difficult to acquire. It is not guaranteed, however, that a good prediction of the shear rheology leads to an accurate description of the extensional stresses too. It might also be possible that the multi-mode Giesekus model employed in this work could not accurately describe the extensional rheology of the PLA. Since at high printing speeds the extensional stresses become more and more important, inaccuracies in the extensional properties would lead to deviations between simulation results and experimental measurements.

4. Conclusions

In this paper, we have investigated with experiments and fluid dynamic simulations the die swell phenomenon observed in a 3D printing process. We have developed an automatic technique to measure the strand diameter in both time and space encoded in a easy-to-use script. Such a technique accounts for possible errors that affect the measurement accuracy due to non-constant strand diameters (induced, for instance, by transient effects) and possible oscillations of the strand with respect to the extrusion direction. A filament of PLA has been extruded at several velocities and three temperatures (with a DoE that includes 48 experiments). We have compared the experimental measurements with numerical results obtained by solving the fluid dynamics equations of the polymer flowing in the die and exiting to the external environment. We have employed an open-source code to simulate the dynamics of the multiphase system, with the PLA shear rheology described by a multi-mode Giesekus constitutive equation.
The experimental analysis method turned out to be accurate and useful to quickly control the extent of the swelling at different velocities and temperatures, opening up the possibility to implement the developed script in an in-line control loop. The simulations describe the die swelling ratio for low and moderate printing velocities well at all three investigated temperatures. At high velocities, however, the experimentally measured swelling ratio suddenly increases and the simulations are not able to capture such a phenomenon. This is likely due to thermal effects that give rise to important temperature inhomogeneities within the filament, not accounted for in the modeling. Furthermore, a better description of the extensional rheological properties might be required, possibly selecting more accurate constitutive equations to describe the PLA rheology. The measurement technique and the numerical method employed in this work can be used to investigate the effect of other relevant parameters on the die swelling phenomenon in 3D printing processes, such as the nozzle geometry and use of different polymers. This will be part of future work.

Author Contributions

Conceptualization, D.T. and G.D.; methodology, G.D. and D.T.; software, D.T. and S.D.R.; validation, D.T. and S.D.R.; investigation, D.T. and S.D.R.; data curation, D.T. and S.D.R.; writing—original draft preparation, S.D.R.; writing—review and editing, G.D. and D.T.; supervision, G.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Turner, B.N.; Strong, R.; Gold, S.A. A review of melt extrusion additive manufacturing processes: I. Process design and modeling. Rapid Prototyp. J. 2014, 20, 192–204. [Google Scholar] [CrossRef]
  2. Gupta, N.; Weber, C.; Newsome, S. Additive Manufacturing: Status and Opportunities; Science and Technology Policy Institute: Washington, DC, USA, 2012. [Google Scholar]
  3. Mackay, M.E. The importance of rheological behavior in the additive manufacturing technique material extrusion. J. Rheol. 2018, 62, 1549–1561. [Google Scholar] [CrossRef]
  4. Serdeczny, M.P.; Comminal, R.; Pedersen, D.B.; Spangenberg, J. Experimental and analytical study of the polymer melt flow through the hotend in material extrusion additive manufacturing. Addit. Manuf. 2020, 32, 100997. [Google Scholar]
  5. Behdani, B.; Senter, M.; Mason, L.; Leu, M.; Park, J. Numerical study on the temperature-dependent viscosity effect on the strand shape in extrusion-based additive manufacturing. J. Manuf. Mater. Process. 2020, 4, 46. [Google Scholar] [CrossRef]
  6. Bourell, D.L.; Leu, M.C.; Rosen, D.W. Roadmap for Additive Manufacturing: Identifying the Future of Freeform Processing; The University of Texas at Austin: Austin, TX, USA, 2009; pp. 11–15. [Google Scholar]
  7. Lombardi, L.; Tammaro, D. Effect of polymer swell in extrusion foaming of low-density polyethylene. Phys. Fluids 2021, 33, 033104. [Google Scholar] [CrossRef]
  8. Polychronopoulos, N.D.; Papathanasiou, T.D. A study on the effect of drawing on extrudate swell in film casting. Appl. Rheol. 2015, 25, 1–7. [Google Scholar] [CrossRef]
  9. Tammaro, D.; Walker, C.; Lombardi, L.; Trommsdorff, U. Effect of extrudate swell on extrusion foam of polyethylene terephthalate. J. Cell. Plast. 2021, 57, 911–925. [Google Scholar] [CrossRef]
  10. Mezi, D.; Ausias, G.; Grohens, Y.; Férec, J. Numerical simulation and modeling of the die swell for fiber suspension flows. J. Non-Newton. Fluid Mech. 2019, 274, 104205. [Google Scholar] [CrossRef]
  11. Robertson, B.; Thompson, R.L.; McLeish, T.C.B.; Robinson, I. Theoretical prediction and experimental measurement of isothermal extrudate swell of monodisperse and bidisperse polystyrenes. J. Rheol. 2017, 61, 931–945. [Google Scholar] [CrossRef] [Green Version]
  12. Tammaro, D. Rheological characterization of complex fluids through a table-top 3D printer. Rheol. Acta 2022, 61, 761–772. [Google Scholar] [CrossRef]
  13. Gifford, W. Compensating for die swell in the design of profile dies. Polym. Eng. Sci. 2003, 43, 1657–1665. [Google Scholar] [CrossRef]
  14. Vlachopoulos, J. Die swell and melt fracture: Effects of molecular weight distribution. Rheol. Acta 1974, 13, 223–227. [Google Scholar] [CrossRef]
  15. Allain, C.; Cloitre, M. Experimental investigation and scaling law analysis of die swell in semi-dilute polymer solutions. J. Non-Newton. Fluid Mech. 1997, 73, 51–66. [Google Scholar] [CrossRef]
  16. Koopmans, R.J. Die swell–molecular structure model for linear polyethylene. J. Polym. Sci. Part A Polym. Chem. 1988, 26, 1157–1164. [Google Scholar] [CrossRef]
  17. Agassant, J.F.; Arda, D.R.; Combeaud, C.; Merten, A.; Muenstedt, H.; Mackley, M.R.; Robert, L.; Vergnes, B. Polymer processing extrusion instabilities and methods for their elimination or minimisation. Int. Polym. Process. 2006, 21, 239–255. [Google Scholar] [CrossRef] [Green Version]
  18. Chang, P.W.; Patten, T.W.; Finlayson, B.A. Collocation and galerkin finite element methods for viscoelastic fluid flow—II. Die swell problems with a free surface. Comput. Fluids 1979, 7, 285–293. [Google Scholar] [CrossRef]
  19. Karagiannis, A.; Hrymak, A.; Vlachopoulos, J. Three-dimensional extrudate swell of creeping Newtonian jets. AIChE J. 1988, 34, 2088–2094. [Google Scholar] [CrossRef]
  20. Bush, M.; Phan-Thien, N. Three dimensional viscous flows with a free surface: Flow out of a long square die. J. Non-Newton. Fluid Mech. 1985, 18, 211–218. [Google Scholar] [CrossRef]
  21. Comminal, R.; Pimenta, F.; Hattel, J.H.; Alves, M.A.; Spangenberg, J. Numerical simulation of the planar extrudate swell of pseudoplastic and viscoelastic fluids with the streamfunction and the VOF methods. J. Non-Newton. Fluid Mech. 2018, 252, 1–18. [Google Scholar] [CrossRef]
  22. Al-Muslimawi, A.; Tamaddon-Jahromi, H.; Webster, M. Simulation of viscoelastic and viscoelastoplastic die-swell flows. J. Non-Newton. Fluid Mech. 2013, 191, 45–56. [Google Scholar] [CrossRef]
  23. Spanjaards, M.; Hulsen, M.; Anderson, P. Transient 3D finite element method for predicting extrudate swell of domains containing sharp edges. J. Non-Newton. Fluid Mech. 2019, 270, 79–95. [Google Scholar] [CrossRef]
  24. Tang, D.; Marchesini, F.H.; Cardon, L.; D’hooge, D.R. Three-dimensional flow simulations for polymer extrudate swell out of slit dies from low to high aspect ratios. Phys. Fluids 2019, 31, 093103. [Google Scholar] [CrossRef]
  25. Spanjaards, M.; Hulsen, M.; Anderson, P. Computational analysis of the extrudate shape of three-dimensional viscoelastic, non-isothermal extrusion flows. J. Non-Newton. Fluid Mech. 2020, 282, 104310. [Google Scholar] [CrossRef]
  26. Spanjaards, M.; Peters, G.; Hulsen, M.; Anderson, P. Numerical study of the effect of thixotropy on extrudate swell. Polymers 2021, 13, 4383. [Google Scholar] [CrossRef] [PubMed]
  27. Azaiez, J.; Guenette, R.; Ait-Kadi, A. Entry flow calculations using multi-mode models. J. Non-Newton. Fluid Mech. 1996, 66, 271–281. [Google Scholar] [CrossRef]
  28. Tang, D.; Marchesini, F.H.; D’hooge, D.R.; Cardon, L. Isothermal flow of neat polypropylene through a slit die and its die swell: Bridging experiments and 3D numerical simulations. J. Non-Newton. Fluid Mech. 2019, 266, 33–45. [Google Scholar] [CrossRef]
  29. Mu, Y.; Zhao, G.; Wu, X.; Zhai, J. Finite-Element Simulation of Polymer Flow and Extrudate Swell Through Hollow Profile Extrusion Die with the Multimode Differential Viscoelastic Model. Adv. Polym. Technol. 2013, 32, E1–E19. [Google Scholar] [CrossRef]
  30. Béraudo, C.; Fortin, A.; Coupez, T.; Demay, Y.; Vergnes, B.; Agassant, J. A finite element method for computing the flow of multi-mode viscoelastic fluids: Comparison with experiments. J. Non-Newton. Fluid Mech. 1998, 75, 1–23. [Google Scholar] [CrossRef]
  31. Mitchell, M.K.; Hirt, D.E. Degradation of PLA fibers at elevated temperature and humidity. Polym. Eng. Sci. 2015, 55, 1652–1660. [Google Scholar] [CrossRef]
  32. Giesekus, H. A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility. J. Non-Newton. Fluid Mech. 1982, 11, 69–109. [Google Scholar] [CrossRef]
  33. Larson, R.G. Constitutive Equations for Polymer Melts and Solutions: Butterworths Series in Chemical Engineering; Butterworth-Heinemann: Waltham, MA, USA, 2013. [Google Scholar]
  34. Hulsen, M.A.; van der Zanden, J. Numerical simulation of contraction flows using a multi-mode Giesekus model. J. Non-Newton. Fluid Mech. 1991, 38, 183–221. [Google Scholar] [CrossRef]
  35. Tsai, W.C.; Miller, G.H. Numerical simulations of viscoelastic flow in complex geometries using a multi-mode Giesekus model. J. Non-Newton. Fluid Mech. 2014, 210, 29–40. [Google Scholar] [CrossRef]
  36. Yao, M.; McKinley, G.H.; Debbaut, B. Extensional deformation, stress relaxation and necking failure of viscoelastic filaments. J. Non-Newton. Fluid Mech. 1998, 79, 469–501. [Google Scholar] [CrossRef] [Green Version]
  37. Pimenta, F.; Alves, M. Stabilization of an open-source finite-volume solver for viscoelastic fluid flows. J. Non-Newton. Fluid Mech. 2017, 239, 85–104. [Google Scholar] [CrossRef]
  38. Pimenta, F.; Alves, M.A. rheoTool User Guide v.6.0. 2022. Available online: https://github.com/fppimenta/rheoTool/blob/master/doc/user_guide.pdf (accessed on 10 October 2022).
  39. Bellini, A. Fused Deposition of Ceramics: A Comprehensive Experimental, Analytical and Computational Study of Material Behavior, Fabrication Process and Equipment Design; Drexel University: Philadelphia, PA, USA, 2002. [Google Scholar]
  40. Go, J.; Schiffres, S.N.; Stevens, A.G.; Hart, A.J. Rate limits of additive manufacturing by fused filament fabrication and guidelines for high-throughput system design. Addit. Manuf. 2017, 16, 1–11. [Google Scholar] [CrossRef]
  41. Osswald, T.A.; Puentes, J.; Kattinger, J. Fused filament fabrication melting model. Addit. Manuf. 2018, 22, 51–59. [Google Scholar] [CrossRef]
  42. Peng, F.; Vogt, B.D.; Cakmak, M. Complex flow and temperature history during melt extrusion in material extrusion additive manufacturing. Addit. Manuf. 2018, 22, 197–206. [Google Scholar] [CrossRef]
Figure 1. Sketch of the 3D printing setup. (a) Feeding zone: setting flow rate; (b) Melting zone: setting temperature; (c) Exit zone: die swell measurement.
Figure 1. Sketch of the 3D printing setup. (a) Feeding zone: setting flow rate; (b) Melting zone: setting temperature; (c) Exit zone: die swell measurement.
Micromachines 14 00329 g001
Figure 2. Geometry of the nozzle used for the experiments.
Figure 2. Geometry of the nozzle used for the experiments.
Micromachines 14 00329 g002
Figure 3. Procedure for die swell measurements. (a) Real image showing a snapshot of the recorder movie. (b) Binarized image with a box selected by the operator. (c) Cropped image from which the horizontal (green dashed line) and perperndicular (yellow dashed line) diameters, D h and D p , are calculated.
Figure 3. Procedure for die swell measurements. (a) Real image showing a snapshot of the recorder movie. (b) Binarized image with a box selected by the operator. (c) Cropped image from which the horizontal (green dashed line) and perperndicular (yellow dashed line) diameters, D h and D p , are calculated.
Micromachines 14 00329 g003
Figure 4. Elastic (red) and viscous (blue) moduli as a function of the frequency (top row) and complex viscosity as a funtion of the shear rate (bottom row) for the three temperatures. The symbols are experimental measurements and the curves are the predictions of the 4-mode Giesekus constitutive equation.
Figure 4. Elastic (red) and viscous (blue) moduli as a function of the frequency (top row) and complex viscosity as a funtion of the shear rate (bottom row) for the three temperatures. The symbols are experimental measurements and the curves are the predictions of the 4-mode Giesekus constitutive equation.
Micromachines 14 00329 g004
Figure 5. Schematic representation of the simulation domain.
Figure 5. Schematic representation of the simulation domain.
Micromachines 14 00329 g005
Figure 6. Initial conditions for the simulations carried out in this work. The domains occupied by the PLA and air are in red and blue, respectively.
Figure 6. Initial conditions for the simulations carried out in this work. The domains occupied by the PLA and air are in red and blue, respectively.
Micromachines 14 00329 g006
Figure 7. Profile of the extrudate for the four different meshes. The z-coordinate is normalized by the inlet diameter, D in . The origin corresponds to the end point of the nozzle.
Figure 7. Profile of the extrudate for the four different meshes. The z-coordinate is normalized by the inlet diameter, D in . The origin corresponds to the end point of the nozzle.
Micromachines 14 00329 g007
Figure 8. (a) Mesh M3 used in the simulations. (b) Zoomed-in image of the die and the exit zone (bottom).
Figure 8. (a) Mesh M3 used in the simulations. (b) Zoomed-in image of the die and the exit zone (bottom).
Micromachines 14 00329 g008
Figure 9. Horizontal and perpendicular strand diameters, D h and D p , measured from the custom-made script, as function of the length of the cropped image: (a) straight strand and (b) bent strand. The two images are taken at different instants of the same experiment at 160 °C and 100 mm/min.
Figure 9. Horizontal and perpendicular strand diameters, D h and D p , measured from the custom-made script, as function of the length of the cropped image: (a) straight strand and (b) bent strand. The two images are taken at different instants of the same experiment at 160 °C and 100 mm/min.
Micromachines 14 00329 g009
Figure 10. Time evolution of the normalized horizontal (blue) and perpendicular (red) diameters for an experiment at 160 °C and 100 mm/min. The bands around the data represent the standard deviation.
Figure 10. Time evolution of the normalized horizontal (blue) and perpendicular (red) diameters for an experiment at 160 °C and 100 mm/min. The bands around the data represent the standard deviation.
Micromachines 14 00329 g010
Figure 11. Axial (left) and radial (right) components of the polymer velocity for three sets of simulation parameters.
Figure 11. Axial (left) and radial (right) components of the polymer velocity for three sets of simulation parameters.
Micromachines 14 00329 g011
Figure 12. z z (left), r r (middle) and r z -component (right) of the stress tensor for three set of simulation parameters.
Figure 12. z z (left), r r (middle) and r z -component (right) of the stress tensor for three set of simulation parameters.
Micromachines 14 00329 g012
Figure 13. Die swell as function of velocity at three temperatures. Comparison between experimental (circles) and simulation results (stars).
Figure 13. Die swell as function of velocity at three temperatures. Comparison between experimental (circles) and simulation results (stars).
Micromachines 14 00329 g013
Table 1. Dimensions (in mm) of the geometrical parameters used in the simulations.
Table 1. Dimensions (in mm) of the geometrical parameters used in the simulations.
L in 3.1 R in 1
L conv 1.2 R die 0.3
L die 1.2 R air 1.1
L air 6 R ext 0.75
Table 2. Values of the constitutive parameters for the 4-mode Giesekus model at the three temperatures.
Table 2. Values of the constitutive parameters for the 4-mode Giesekus model at the three temperatures.
T = 160 °C λ i (s) α i η i (Pa · s)
Mode 10.2640.5811,556
Mode 23.130.445244.4
Mode 30.03080.516943.7
Mode 40.0028230.501122.9
T = 180 °C λ i (s) α i η i (Pa·s)
Mode 10.02250.503124
Mode 20.1910.672959
Mode 32.4070.14791.4
Mode 40.00230.50860.2
T = 200 °C λ i (s) α i η i (Pa·s)
Mode 10.01560.521394
Mode 20.1260.68907.6
Mode 31.5470.10149.8
Mode 40.001730.50578.2
Table 3. Physical parameters for PLA and air.
Table 3. Physical parameters for PLA and air.
ρ γ ρ air η air
(kg/m 3 )(mN/m)(kg/m 3 )(Pa·s)
1000421.2251.8 · 10 5
Table 4. Dimensionless numbers for the three temperatures considering the maximum inlet average velocity investigated in this work.
Table 4. Dimensionless numbers for the three temperatures considering the maximum inlet average velocity investigated in this work.
TemperatureReWiCa
160 °C4.01 · 10 7 7.802.96 · 10 3
180 °C1.29 · 10 6 6.029.21 · 10 3
200 °C3.30 · 10 6 3.873.61 · 10 3
Table 5. Details of the meshes used in the spatial convergence study.
Table 5. Details of the meshes used in the spatial convergence study.
Mesh# of Cells# of Elements on R die
M114286
M2571212
M322,84824
M435,70030
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

De Rosa, S.; Tammaro, D.; D’Avino, G. Experimental and Numerical Investigation of the Die Swell in 3D Printing Processes. Micromachines 2023, 14, 329. https://doi.org/10.3390/mi14020329

AMA Style

De Rosa S, Tammaro D, D’Avino G. Experimental and Numerical Investigation of the Die Swell in 3D Printing Processes. Micromachines. 2023; 14(2):329. https://doi.org/10.3390/mi14020329

Chicago/Turabian Style

De Rosa, Stefano, Daniele Tammaro, and Gaetano D’Avino. 2023. "Experimental and Numerical Investigation of the Die Swell in 3D Printing Processes" Micromachines 14, no. 2: 329. https://doi.org/10.3390/mi14020329

APA Style

De Rosa, S., Tammaro, D., & D’Avino, G. (2023). Experimental and Numerical Investigation of the Die Swell in 3D Printing Processes. Micromachines, 14(2), 329. https://doi.org/10.3390/mi14020329

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop