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.