Fibre Bridging and Nozzle Clogging in 3D Printing of Discontinuous Carbon Fibre Reinforced Polymer Composites: Coupled CFD-DEM Modelling

: A coupled multiphase model based on computational fluid dynamics (CFD) and discrete element method (DEM) is developed to numerically investigate the extrusion-based 3D printing process of discontinuous carbon fibre reinforced polymer composites. Short carbon fibres are modelled as rigid bodies by clumping discrete spheres in DEM, while polymer matrix is treated as an incompressible Newtonian fluid in CFD. A fluid-particle interaction model is adopted to couple DEM and CFD and represent the dynamic fibre/matrix interaction. Collisions between fibres are considered naturally in DEM by using the Hertz-Mindlin contact law. The coupled CFD-DEM is validated, both qualitatively and quantitatively, against X- ray microtomography (μCT) experimental results for the T300/PA6 composite. Parametric study on various fibre lengths, fibre volume fraction and resin viscosity using the CFD-DEM model shows that the nozzle clogging tends to occur when the fibre length and/or the fibre volume fraction are increased. Use of a polymer matrix with a lower viscosity can be effective to eliminate the clogging issue when printing composites with relatively short fibres. The fibre length is dominating when long fibres are used and the clogging is largely independent on the viscosity of the polymer matrix. Finally, a potential solution of using a cone sleeve insert located above the shrinking region to address the nozzle clogging issue is proposed and numerically assessed.


INTRODUCTION
Carbon fibre reinforced polymer (CFRP) composites offer various advantages such as high strength and high stiffness-to-weight ratio, as well as excellent fatigue performance and corrosion resistance [1,2]. Apart from the traditional manufacturing methods, the recently-emerged additive manufacturing techniques (AM, also known as 3D printing) have shown the potential to enable low-cost and mould-free manufacturing of three-dimensional objects with complex geometries [3]. Among them, the extrusion-based fused filament fabrication (FFF) printing is the most popular process (by number of machines), which fabricates a part by melting and depositing the composite material layer by layer [4]. Due to the favourable melt processability, a wide variety of thermoplastic can be used as the matrix materials in the FFF 3D printing of composites, such as ABS, PLA and nylon. Discontinuous carbon fibres are usually introduced to reinforce the thermoplastic with a purpose of achieving printed composites with improved strength, stiffness and heat resistance [5].
Previous researches have shown that improved mechanical properties can be achieved in the printed composites [6], in particular the flexural properties could be improved by 2~3 times with the mixed discontinuous carbon fibres [7].
Thermoplastic polymers have a relatively higher melt viscosity than thermosetting polymers and it is challenging to form strong bonds between a thermoplastic polymer and carbon fibres. The chopped carbon fibres mixed in the melt thermoplastic during FFF printing would further increase the viscosity, and the difference in their coefficients of thermal expansion would result in more voids during the heat cycling [6,8,9]. Although manufacturing of short fibre composites at a faster production rate is desirable with the FFF technique, if not well treated, the mechanical properties of finished parts could be impaired and even become inferior to that of the printed neat polymer materials, which is mainly caused by the voids formed in the printing process [10]. On the other hand, the chopped carbon fibres currently used in the 3D printing are usually shorter than 0.2 mm, which is not long enough to make a distinctive enhancement on the mechanical properties of the printed composites. Markforged ® reported that the tensile properties of discontinuous fibre reinforced thermoplastic composites are lower than the 3D printed neat thermoplastic.
Attempting to address the afore-mentioned issues, the mechanisms of the melt flow as well as void formation through the standard FFF printing process have been investigated [11]. Previous researches have studied the effect of printing parameters on the mechanical performance [12][13][14][15]. And the in-situ consolidation techniques such as ultrasonic vibration [16] and post printing treatment such as hot-pressing [17] have been adopted to remove the air voids induced during the printing. In order to achieve higher mechanical performance, other attempts have been made such as the use of longer discontinuous fibres (e.g. 3 mm) in previous work [18]. A few studies also aimed to customise the alignment of discontinuous fibres in the composites with the nozzle design and paths optimization [19][20][21], in addition to rotational printing [22] and the use of a magnetic field [23]. Another extrusion-based printing technique, direct-ink writing (DIW) has been used to print short fibre reinforced polymer composites, in which a thermosetting resin with low viscosity is mixed with the short carbon fibres prior to the printing [20]. In the DIW printing process, a long and needle-like nozzle is usually adopted, but risks of fibre bridging and nozzle clogging still exist in the regions close to the nozzle outlet when long fibres are used [24].
Due to the lack of understanding on fibre motion and distribution as well as the interaction among them, the mechanism for the issues of fibre breakage and nozzle clogging is still unclear. In addition, the large variability of printing parameters results in inconsistent mechanical properties of the printed composites. Optimal printing parameters are usually identified through a nontrivial trial-and-error process for each specific composite, followed by traditional tensile/flexural tests and the characterisation of fracture profile [12,25]. Therefore, there is still a lack of established process-properties relationship for 3D printing of composites in general.
Computational modelling, as an alternative research tool, can provide much more details of multiphase flow which are hardly accessible via traditional physical measurement and thus has been adopted to improve the understanding of the 3D printing process [26].
Papon et al. [27] developed a computational model and numerically investigated the effects of nozzle speed, feed rate and nozzle geometry on the mechanical performance of 3D printed composites. Finite element method (FEM) was also used to evaluate the polymer melt flow and the associated effects of nozzle geometry and extrudate swell on fibre orientation [28].
The effective bulk mechanical properties, specifically the longitudinal and transverse moduli, and the coefficient of thermal expansion, were calculated for the deposited bead based on the spatially varying fibre orientation tensors [29]. Yang et al., [30] established a new coupling model based on smoothed particle hydrodynamics (SPH) and discrete element method (DEM) to simulate the 3D printing process of short CFRP composites and obtained promising results for the rheological flow and fibre orientation.
From the literature survey, it is still necessary to further understand (and eventually better control) the flow characteristics of the melt composites during the 3D printing process.
Numerical simulations with rigorous validations would be useful to achieve this target.
Unfortunately, modelling such multiphysics, multiphase flow process is complicated and there has been very few studies yet. Croom et al. [24] used computational fluid dynamics (CFD) to simulate the resin flow in DIW 3D printing process. However, fibres were not explicitly considered in the model and fibre-fibre collisions could not be fully represented, thus the nozzle clogging was not well addressed. The fibre-fibre collision behaviours were captured in the SPH-DEM model by Yang et al. [30], but the fibre length is very short and the occurrence of nozzle clogging was investigated due to the fact that the model was two dimensional. Therefore, this study aims to develop a 3D coupled model based on CFD and DEM to understand the effects of printing parameters, including fibre volume fraction (Vf), fibre length and matrix viscosity, on the fibre misalignment and nozzle clogging in extrusion-based 3D printing of composites. The paper is organised as follows. Section 2 introduces the numerical methods for the coupled CFD-DEM model and its validation, followed by comprehensive parametric studies in Section 3. Potential solutions are proposed and numerically assessed in Section 4. Conclusions are finally drawn in Section 5.

CFD-DEM coupling
In the coupled CFD-DEM model, short fibres are regarded as a discrete, solid phase and simulated using DEM, while the polymer matrix is simulated using CFD. As the short fibres normally have an aspect ratio larger than one, each fibre is represented by a group of discrete DEM particles that are bonded together as a rigid clump. The motion of each particle in DEM is described by Newton's second law of motion. For particle i, the equations governing the translational and rotational motions can be expressed as: where mi and Ii are the particle mass and the moment of inertia, respectively. vi and ωi are the linear and angular velocities, respectively. ni is the number of particles interacting with particle i. Relevant forces are gravitational force mig, fluid-particle interaction force Fpf,i and inter-particle collision forces including elastic force Fe,ij and viscous damping force Fd,ij.
The Hertz-Mindlin contact model is adopted to calculate the inter-particle collision force.
The torques acting on particle i by particle j include Tt,ij generated by the tangential force, Tr,ij known as the rolling friction torque and Tn,ij generated when the normal force does not pass through the centre of particle i.
The polymer matrix is regarded as a continuous phase whose motion follows the Navier-Stoke equations given by: where ρf is fluid density, t is time, uf is fluid velocity and p is fluid pressure.
is the voidage where kv is the total number of particles in the current computational cell of volume ΔV and Vi is the volume of particle i. τ is the stress tensor which is given as: The drag force is the only fluid-particle interaction force considered in this study, which can be described as: where β is the fluid-particle interphase drag coefficient [31] and can be defined as: where dp,i is the diameter of particle i. CD is the drag coefficient, which can determine the effective transfer of the interphase force. The formula has considered the effect of the particle shape using the concept of sphericity φ, as given by: 23 where uv is the Reynolds number for particle i. dv is the equivalent diameter defined as the diameter of a sphere with the same volume as the non-spherical particle. Therefore, dp,i = dv in this study.

3D printing experiments and X-ray µCT scans
To validate the computational model, 3D printing experiments of discontinuous carbon fibre reinforced polymer were carried out on a Prusa i3 MK3S+ printer using a printing material Onyx ® from Markforged ® . The material data sheet has shown that Onyx ® filament has diameter of 1.75 mm and consists of polycaprolactam (PA6) and short T300 carbon fibre [32]. The thermal properties of Markforged® ONYX were characterized in the previous study [33]. The glass transition temperatures (Tg), crystallisation temperature (Tc) and melting point (Tm) are 27•C, 161•C and 197•C respectively. In our previous study [34], a X-ray microtomography (µCT) scan of the intact filament was conducted, with an accelerating voltage of 80 kV and a power of 7 W. The exposure time and effective pixel size for each scan can be found in [34]. It demonstrated that it has a mean fibre length of 0.12 mm, an average fibre diameter 7 µm and a fibre volume fraction of about 13.34%. The µCT scan also shows that the short carbon fibres are uniformly distributed and most of them are aligned with the longitudinal direction of the filament. The tensor component aZZ is about 0.9 in the raw filament (z-axis direction is the flow direction). In the printing experiment, a temperature of 275°C was used as suggested by Markforged ® , together with a filament feed rate of 0.0012 m/s.
The internal geometry of the nozzle is essential for the configuration of the computational model, so a µCT scan of the nozzle was also carried out, as shown in Fig.   1a. It can be seen the inner diameter of the nozzle is 2 mm, and it is reduced to 0.45 mm in the outlet with a certain slope, as indicated in Fig. 1b. A composite sample was taken from the nozzle (see Fig. 1c) and scanned to obtain the microscopic information of fiber distributions, which will be used for the validation of the computational model. More details of the printing experiments can be referred to our previous work [34].  The main parameters used in the numerical simulations are listed in Table 1. Identical to the experimental settings, the viscosity of the PA6 resin is set as μf=120 Pas, and the vertically aligned fibres are gradually and randomly generated until reaching a volume fraction Vf of 13.34%. areas [34]. The local excess of matrix in section Ⅰ may result from both shear re-orientation between core and edge and the strong re-orientation at the entrance [35]. It also implies that the fibre fraction of the central region is initially close to the nominal value but    Corresponding numerical results are summarised in Table 2. If the simulation is successfully completed with continual flow of fibres through the nozzle outlet, it is considered that there is no fibre bridging and thus no nozzle clogging. On the contrary, if the fibres are gradually clogging in the nozzle and only polymer resin flows through the nozzle outlet, the simulation is terminated (due to an abrupt change of the concentration of DEM particles in the sharply narrowed region) and it is considered as an occurrence of fibre bridging and nozzle clogging.
As can be seen in length and fibre volume fraction, and there would be a scope for optimising the combination of these two factors to mitigate the risk of nozzle clogging. Therefore, a more detailed and quantitative study of the fibre flow through the entire printing process is carried out in the next section.   A typical bridging phenomenon in the evolution of fibre flow for L=0.24 mm and initial Vf=26.68% is given in Fig. 6. It can be seen that as the fibres enter the shrinking region, the fibre orientation in area III is changed substantially from vertical to inclining.

Effects of fibre length and fibre volume fraction
Severe interlocking takes place before these inclining fibres entering area V. Once the interlocking structure is strong enough to support the weight of above fibres, fibre bridging is formed and no more fibres can pass through the nozzle aperture to the nozzle outlet. As a result, the 3D printing process is terminated due to the clogging in the nozzle and the sudden rise of pressure in nozzle inlet, which prevents further feeding of the filament.
Therefore, the current numerical results recommend that it is critical to control the fibre orientation in the shrinking region to prevent fibre bridging. The occurrence of bridging is related to the average fibre (i.e. DEM particle) velocity in the nozzle as shown in Fig. 7. Recall that each fibre is formed by a group of spheres thus the average particle velocity is a statistical representation of the fibre flow. It can be observed from Fig. 7a that average particle velocities in all the cases show quite similar trends, i.e. increasing at first (due to gravity), then dropping (due to the obstruction of the shrinking section) and finally reaching an almost constant value (due to the narrow section in the nozzle outlet). In this study, nozzle clogging is regarded to take place when the constant value is lower than 0.5 mm/s. It is seen that the 3D printing processes can run well when L=0.12 mm regardless of the initial Vf. This means that adopting shorter fibres is beneficial for the fluidity at the same initial Vf, and would not cause fibre bridging even with high initial Vf. When fibre length L=0.24 mm is used as shown in Fig. 7b, it can be seen that the particle flow is not terminated with the increase of the initial Vf until it reaches 26.68% and the average particle velocity becomes very low afterwards. Low velocities are commonly observed when the fibre length further increases to L=0.35 mm as shown in Fig. 7c. Long fibres flow very slow even with a low initial Vf. Therefore, it can be concluded that the fibre length and the initial Vf have a combined influence on the 3D printing process. Since a long fibre and high initial Vf are favored, the current numerical results suggest to adopt the case of L=0.24 mm with initial Vf=20.01% or L=0. 35  bridging actually takes place in areas Ⅲ and Ⅳ, so it is necessary to carry out a microscale analysis on the particle structure in these areas to further investigate the mechanisms. For the same reason, the analysis is only carried out for Markforged PA6 but all the findings would be valid for Arkema Elium 150 and Victrex PEEK 381G. Vf=26.68%. The particle coloured by their CN. The cylinders between particles represent the particle-particle contacts.

Microscale analysis on the particle structure
Based on the observations in Fig. 6, area Ⅳ is of interest and the distribution of Coordination Number (CN) in this area is shown in Fig. 8. CN of a DEM particle is defined as the number of surrounding particles contacting with it. Therefore, overall distribution of CN is a good indicator to show the packing structure of particles. Fig. 8a presents the probability density function (PDF) of CN in area Ⅳ for different study cases, from which it can be seen that the CN curves for all these cases show a Gaussian-like profile. For each fibre length (and therefore different number of spheres for the fibre), with the increase of the initial Vf, the CN curves shift parallel to the larger CN slightly. This is reasonable because an increase of the initial Vf leads to more potential contacts in a constrained space. For a given initial Vf, it is interesting to see that the peak positions of these CN curves are not influenced by the increase of the fibre length. Normally, longer fibres would have larger CNs in a dense packing system due to their larger surfaces for contacting with more other particles. However, that seems not the actual case here mainly  When the contact structure of the packed particles is obtained, it is useful to find out the 21 contact forces at each contact point. Fig. 9 shows the PDF of the contact forces and their fittings. It can be seen that the overall distribution trends of the contact forces are quite similar for all considered initial Vf and fibre length. That is, most of the contact forces (i.e. >60%) are comparable and only a small proportion of larger forces exist. However, it is clear that the magnitude of the contact forces increases with both initial Vf and fibre length. This is because the mass of a single fibre increases with its length. Larger particle mass leads to larger contact forces. Furthermore, there are more fibres with the increase of the initial Vf and thus the packing structure become much denser, resulting larger contact forces. Above microscale analysis shows that the fibre distribution in area III is more sensitive to the fibre length and initial Vf. Therefore, to prevent the formation of fibre bridging and nozzle clogging, a cone sleeve insert is installed in area III, which is located in the nozzle aperture right above the shrinking region, as shown in Fig. 10.  Table 3. It can be seen that nozzle clogging is successfully resolved by Plan B rather than Plan A. The average particle velocities for both designs are shown in Fig. 11 for different initial Vf. As shown, for both the two designs with two initial Vf, the general trends of the average particle velocity curves are quite similar with the original one. The main discrepancy takes place at the 'dropping' state when the fibres enter the nozzle aperture. Due to the existence of the insert, the average particle velocity drops faster than that in the original nozzle. This is especially true for Plan A, in which the insert is found to retard the fibre sedimentation greatly and nozzle clogging cannot be avoided.

Proposed solutions to avoid nozzle clogging
Therefore, only when the cone sleeve insert is properly designed such as Plan B, the long fibres can flow into area V without any bridging taking place unlike the original nozzle and Plan A.

Concluding remarks
In this work, the extrusion-based 3D printing process of carbon fibre reinforced polymer composites with three different fibre lengths and three different polymer resins considered is numerically investigated using a coupled CFD-DEM approach. Effects of fibre length, fibre volume fraction and polymer property on the 3D printing process are discussed and the occurrence of nozzle clogging is predicted. Microscale analysis on the fibre particle structure is conducted. Finally, potential solutions to avoid nozzle clogging are proposed and numerically evaluated. According to the numerical results, following conclusions can be drawn for the range of parameters considered in the current study.
(1) Adopting short fibre is profitable for the fluidity and avoiding nozzle clogging at the same initial Vf. At a relatively short fibre length L=0. 24  Greek letters     Comparison of bre volume fractions in different representative areas.  Evolution of bre ow for L=0.24 mm and initial Vf=26.68% (only bres in the shrinking region are shown).

Figure 7
The average particle velocity in the system with different bre length and initial bre volume fraction when μf=120 Pa.s and ρ=1150 kg/m3.  Distribution of CN in area for Markforged PA6: (a) tting lines of the PDF of CN, (b) distribution of CN for Vf=13.34%, (c) distribution of CN for Vf=20.01% and (d) distribution of CN for Vf=26.68%. The particle coloured by their CN. The cylinders between particles represent the particle-particle contacts.  Sketches of different cone sleeve inserts in the nozzle.

Figure 12
Distributions of DEM particles (top) and DEM particles colored by their contact forces (bottom) for L=0.35 mm and initial Vf=26.68% in different nozzles.

Figure 13
The percentage of the orientations of bre particles for L=0.35 mm and initial Vf=26.68% in different nozzles.