2D Fusion Simulations and Experimental Confirmations of Print Paths Using Composite Particles with Particle Method for Fused Filament Fabrication

Yuto Imaeda¹, Akira Todoroki¹, Ryosuke Matsuzaki², Masahito Ueda³, Yoshiyasu Hirano⁴

¹Department of Mechanical Engineering, Tokyo Institute of Technology, Meguro, Japan.
²Department of Mechanical Engineering, Tokyo University of Science, Noda, Japan.
³Department of Mechanical Engineering, Nihon University, Chiyoda, Japan.
⁴Japan Aerospace Exploration Agency, Mitaka, Japan.

DOI: 10.4236/ojcm.2022.124009                

Abstract

Printing short fibre/thermoplastic composites using the fused filament fabrication method sometimes creates a gap between print paths. In this study, the two-dimensional moving particle semi-implicit method for liquid simulation was applied to simulate the print-path fusion process. The three-dimensional movement of the nozzle was simulated using the sliding motion of the nozzle. The method was applied to the printing of short carbon fibre/polyamide-6 composites, and the simulation results were compared with those of experiments. The simulated results of the cross-sectional configuration agreed well with the experimental results. This will enable the optimization of printing process parameters thus reducing the gap between print paths.

Keywords

Carbon Fiber, Thermoplastic Resin, Computational Modelling, Short Carbon Fiber, Three-Dimensional Printer, Fusion Process

1. Introduction

The fused filament fabrication (FFF) method is widely adopted for thermoplastic three-dimensional (3D) printing. Thermoplastic FFF-type 3D printers have been mainly used for prototyping in the engineering fields. In this technology, thermoplastic filaments with short fibers are generally employed as filaments owing to their high stiffness and strength.

Love et al. [1] showed the high stiffness and dimensional stability of the 3D printed products using short carbon fibre/acrylonitrile-butadiene-styrene composite filaments. Tekinalp [2] et al. revealed the high strength of the short carbon fibre/ABS composites; the fibers in the printed path were shown to be aligned in the printing direction. Postiglione et al. [3] added carbon nanotubes to polylactic acid (PLA) and investigated the electrical conductivity of the resulting composites. They showed that electrically conductive plastic products could be fabricated using the FFF-type 3D printers. Matsuzaki [4] et al. fabricated continuous carbon-fibre-reinforced PLA. The use of continuous carbon fibers enabled the printing of high-strength and high-stiffness products using the FFF-type 3D printers. Wang et al. [5] and Parandoush and Lin [6] reviewed research papers on 3D printed composites.

Tekinalp et al. [2] conducted tensile tests with various fibre weight fractions using a short carbon fiber/ABS composite. With an increase in the number of carbon fibers, the number of voids increased in the cross-section of the specimens. Ning et al. [7] performed tensile tests on short carbon fiber/ABS composites and found an increase in voids with an increase in the fiber volume fraction. Ferreira et al. [8] performed tensile tests on short carbon fiber/PLA composites and demonstrated fibre alignment in the print path orientation. Ning [9] investigated the effect of the process parameters on the tensile properties of short carbon fibre/ABS composites. They found that the optimal nozzle temperature and infill speed maximized the tensile strength. Zhang et al. [10] conducted tensile and shear tests using a short carbon fibre/ABS composite. They found the largest gap between the print paths in the ±45˚ specimens. Naranjo-Lozada [11] investigated the effect of the infill pattern of the 3D printer, Mark Two. Yavas et al. [12] performed tensile and fracture toughness tests on a 3D printed cross-layer short carbon fiber/polyamide 6 (PA-6) Onyx. The tensile test results showed that the strength of the layup direction was approximately 32% of the in-plane strength. However, the effects of bed temperature and extrusion rate were not reported [13]. These studies indicated that the 3D printed composites require optimization of the fabrication conditions to reduce the gap between the print paths or voids. However, experimental optimizations are time-consuming, and an appropriate fabrication process simulation method is required.

Talagani et al. [14] used the finite element method (FEM) to calculate the residual stress and predict cracking in the printing of a full-size working car. They used short carbon/ABS composites, and cracking was successfully predicted. Makino et al. [15] adopted the smoothed particle hydrodynamics (SPH) method to compute resin flow during the printing process of the FFF method. They calculated the fusion process of thermoplastic resin using thermal conduction analysis, although they did not adopt latent heat for the phase change. Yang et al. [16] performed an SPH simulation of the 3D printing of short carbon fibre/PA-6 composites. The short carbon fibers were modelled using a discrete element method. The study dealt with the extrusion of PA-6 particles with carbon fibers from a nozzle at a constant temperature. Bertevas et al. [17] computed the shear deformation at the nozzle using the SPH method. They used rods made from rigid particles as fibers at constant temperature. The fibre direction during extrusion from the nozzle was dealt with in the study. Brenken et al. [18] calculated the heat transfer and shrinkage of the fabrication process during 3D printing of short carbon fibre/polyphenylene sulfide (PPS) composites. They used FEM analysis, and the anisotropic material properties were considered. The crystallization of the PPS resin was also calculated in the study. Ouyang et al. [19] used the SPH method to investigate the effect of the Péclet number of the heat transfer on the nozzle. Rigid rods were used to represent short fibers. Although the latent heat from liquid to solid was not considered in this analysis, the heat transfer and fibre alignments in the nozzle were carefully investigated. Ouyang et al. [20] also considered the effect of the anisotropic thermal conductivity caused by fibre alignments in the SPH analysis. Zhang et al. [21] calculated the nozzle coggling using computer fluid dynamics and discrete element method. Lee etal. [22] used finite element method to compute the thermoplastic flow during printing process. Farazin and Mohammadimehr [23] used molecular dynamics to simulate the 3D printing process.

Several researchers [16] - [23] have simulated the fabrication process of 3D printing of short fibre/thermoplastic composites. They used the FEM analysis or the SPH analysis method. FEM analysis requires the target structure to be divided into meshes. Thus, thermoplastic resin flow analysis is difficult to perform because it requires remeshing of the structure. SPH analysis uses particles as nodes of an entire analysis field. The movements and properties of the target field are expressed by the response surface of the radial basis functions (kernel functions) at the nodes [24]. Thus, the SPH method does not require remeshing to analyze the resin flow. In studies that have used the SPH method, the fibers are modelled as serially linked particles. Thus, when we adopt an actual sized nozzle and fibers, a large number of particles are required, which results in considerable computational time. Moreover, many studies that have used the SPH method adopt an explicit analysis method. This requires a small-time step when the particle size is small to satisfy the Courant-Friedrichs-Lewy condition. A small-time step generally indicates a long computational time.

We adopted the moving particle semi-implicit (MPS) method for thermoplastic flow analysis during 3D printing [25]. In the MPS method, the particle is considered to have spherical properties; a gradient model and a Laplacian model were proposed to calculate fluid from the discrete values around the particle, and a semi-implicit method was used to calculate the incompressible flow [25].

In our previous study, a composite particle was adopted, and the solidification process of the printed path was calculated by considering the latent heat and non-linear viscosity of the resin; the results were compared with experimental results [25]. In the previous study, the MPS method was applied using a 3D model that required a large number of particles, although it was a single-path analysis. Therefore, it was difficult to perform simulations of the fusion processes between the printed paths that required more particles. In the present study, a new 2D analysis is proposed to simulate 3D nozzle movements. First, single-path printing was simulated to confirm the newly developed 2D simulations of the 3D printing nozzle movements. Second, the method was applied to the print path fusion process, and the results were compared with those of experiments

2. MPS Method

2.1. Modelling of Melting and Solidification Using Composite Particles

The melting and solidification of thermoplastic resin with short fibers were modelled to simulate the heat transfer and flow phenomena during the 3D printing process in a previous study [26]. This method is briefly explained in this section.

Voller et al. [27] developed a thermofluidic simulation method to model the melting and/or solidification process with flow approximation in porous solids. In this approximation, the porous solid with melted thermoplastic has porosity and enthalpy. The porosity η is defined as a value that linearly increases with the phase change from the beginning to the end of the melting process. It can be expressed as follows:

$\eta =\{\begin{array}{ll}0\hfill & \left(H<H_S\right)\hfill \\ \frac{H-H_S}{H_S-H_L}\hfill & \left(H_S\le H\le H_L\right)\hfill \\ 0\hfill & \left(H_S<H\right)\hfill \end{array}$ (1)

where H_s is the enthalpy at the beginning of melting and H_L is the enthalpy at the end of melting. The temperature T is defined as follows:

$T=\{\begin{array}{ll}{T}_S+\frac{H-H_S}{\rho C_p}\hfill & \left(H<H_S\right)\hfill \\ T_S\hfill & \left(H_S\le H\le H_L\right)\hfill \\ T_S+\frac{H-H_L}{\rho C_p}\hfill & \left(H_L<H\right)\hfill \end{array}$ (2)

where T_s is the melting point and C_p is the specific heat after complete melting.

Enthalpy H is the sum of the sensible heat h and latent heat ΔH. The latent heat is released linearly with the increase in porosity as follows:

$H=h+\Delta H$ (3)

$\Delta H={\eta }_{\text{L}}\left(H_L-H_S\right)$ (4)

This method approximates the material with the liquid and solid phases during melting or solidifying as a flow in a porous medium. This approximation method is adopted in the simulation of the phase change of the metal. In Equation (4), η_L is the liquid fraction. When the enthalpy equation, Equation (4), is applied to the composite materials, we should consider the extended liquid fraction η_c for composite materials, which includes the fibre volume fraction. The conventional definition of the liquid fraction of metal η_m, which is the volume fraction of the liquid phase, is as follows:

${\eta }_m=1-F_s$ (5)

where F_s is the solid-phase fraction of the metal. To express a composite particle that comprises fibers, liquid thermoplastic, and solid plastic, the definition of the porous rate is modified as follows:

${\eta }_c=1-\left(V_f+F_s\right)$ (6)

where η_c is the liquid fraction of the particle, V_f is the fiber volume fraction, and F_s is the solid-phase fraction of the thermoplastic. The relationship between temperature, liquid fraction, and enthalpy is shown in Figure 1.

In Figure 1, the abscissa is the enthalpy, the left ordinate is the liquid fraction, and the right ordinate is the temperature. When the enthalpy is lower than the value when melting starts (H_solid), the composite particle comprises solid thermoplastic and fibers. The temperature increases linearly with the increase in heat flow. When the enthalpy is higher than H_solid and lower than the perfect liquid value of thermoplastic (H_liquid), the temperature is constant, and the heat flow is used to transform solid thermoplastic to liquid thermoplastic. When the enthalpy is higher than H_liquid, the thermoplastic is liquid, and the heat flow is used to increase the temperature of the liquid thermoplastic and fibers.

2.2. Simulation Method with MPS

In the present study, the MPS method [26], which was developed for the implicit analysis of incompressible fluid flows, was adopted for the short fibre composite flow simulation of the FFF process. In the MPS method, remeshing of the FFF process of a short fibre/thermoplastic composite flow is not required. Although carbon fibers are highly oriented in the printing direction [2], the present study regarded short carbon fibre thermoplastic composites as isotropic materials because 2D analysis perpendicular to the print path was considered in this study.

The material properties of the particle in this study under the MPS method were obtained from those of the composite materials comprising short fibers and thermoplastic. The material properties of the composites can be calculated using the rule of mixtures. For example, density can be obtained as follows:

${\rho }_c={\rho }_f{V}_f+{\rho }_m{V}_m$ (7)

where ρ_c, ρ_f, and ρ_m are the densities of the composite, fiber, and matrix, respectively. The subscripts f and m denote the material property of the fiber and matrix, respectively. Similarly, the subscript c indicates the material properties of the composites. In the MPS method, the fluid flow can be expressed using the equation of continuous and the equation of motion as follows:

$\frac{D\rho }{Dt}+\rho \nabla \cdot u=0$ (8)

$\frac{Du}{Dt}=-\frac{1}{\rho }\nabla P+\nu {\nabla }^{2}u+g$ (9)

Figure 1. Relationship between liquid fraction, temperature, and enthalpy in a composite particle model with phase change.

where D/Dt is the Lagrange derivative, u is the velocity vector, ρ is the density, ν is the kinematic viscosity coefficient, and g is the vector of gravitational acceleration. Using the interparticle interaction model shown in Figure 2, the derivative and Laplacian in Equation (9) can be discretized [26].

In discretization, a weight function is used to express the interaction between particles #i and #j, where particle #j is located at a distance of r from particle #i.

$w\left(r\right)=\{\begin{array}{ll}\frac{r_e}{r}-1\hfill & \left(r<r_e\right)\hfill \\ 0\hfill & \left(r\ge r_e\right)\hfill \end{array}$ (10)

where w(r) is the weight function and r_e is the effective radius of the interaction zone. The density in Eqs. (8) and (9) is replaced by the density of particles n_i as follows:

$n_i={\displaystyle \underset{j\ne i}{\sum }w\left(\left|r_j-r_i\right|\right)}$ (11)

where $w\left(\left|r_j-r_i\right|\right)$ is the weight function. Using the interparticle interaction mode, the discrete model of the gradient can be expressed as follows:

${\langle \nabla \varphi \rangle }_i=\frac{d}{n^0}{\displaystyle \underset{j\ne i}{\sum }\frac{{\varphi }_j-{\varphi }_i}{{\left|r_j-r_i\right|}^2}}\left(r_j-r_i\right)w\left(\left|r_j-r_i\right|\right)$ (12)

where ϕ_i is the physical quantity, such as pressure, temperature, etc. of particle #i, d is the dimension of the target space, and n⁰ is the density of particles at time t = 0. The discrete model of the Laplacian is as follows:

${\langle {\nabla }^2\varphi \rangle }_i=\frac{2d}{{\lambda }^0{n}^0}{\displaystyle \underset{j\ne i}{\sum }\left({\varphi }_j-{\varphi }_i\right)w\left(\left|r_j-r_i\right|\right)}$ (13)

${\lambda }^0=\frac{{\displaystyle \underset{j\ne i^{\prime }}{\sum }{\left|r_j^0-r_i^0\right|}^{2}w\left(\Vert r_j^0-r_i^0\Vert \right)}}{{\displaystyle \underset{j\ne i^{\prime }}{\sum }w\left(\Vert r_j^0-r_i^0\Vert \right)}}$ (14)

Figure 2. Particle interaction model.

where λ⁰ is the normalization constant of the Laplacian model, and r⁰ is the initial equal spacing between the particles. The explicit calculation of the viscosity and gravity terms in Equation (9) gives the tentative velocity of the target particle. The tentative velocity is modified after the calculation of pressure. In this study, pressure can be calculated using Tanaka’s method [28] implicitly as follows:

${\langle {\nabla }^{2}P\rangle }_i^{k+1}=\frac{{\rho }_0}{\Delta t}{\langle \nabla \cdot u\rangle }_i^{*}+\delta \frac{{\rho }_0}{{\left(\Delta t\right)}^2}\left(\frac{n_i^{*}-n^0}{n^0}\right)$ (15)

where $n_i^{*}$ is the density of particle #i at the tentative location r^*, ρ₀ is the initial density, Δt is the time step, superscript k denotes the time step in the numerical computations, and δ is a computational parameter usually set to a positive real number smaller than 1. u^* is the tentative velocity, which is calculated as follows:

$u_i^{*}=u_i^k+\left[\frac{2d}{n^0}{\nu }_{ij}{\displaystyle \underset{j\ne i}{\sum }\frac{\left(u_j-u_i\right)}{{\left|r_j-r_i\right|}^2}}\left(r_j-r_i\right)w\left(\left|r_j-r_i\right|\right)+g\right]\Delta t$ (16)

The calculated pressure is substituted into Equation (9). This gives the true velocity of the next time step.

The following heat energy conservation law is used in the heat transfer calculation in the MPS method.

$\frac{DH}{Dt}=\lambda {\nabla }^{2}T$ (17)

where H is enthalpy, λ is the thermal conductivity, and T is the temperature.

The viscosity of a short fiber thermoplastic is represented by the Carreau model [29] as follows:

$\mu ={\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left\{1+{\left(\beta \stackrel{\dot{}}{\gamma }\right)}^2\right\}}^{\frac{n-1}{2}}$ (18)

where μ₀ and μ_∞ are zero and infinite shear rate viscosities, respectively, β is the characteristic time, n is the material constant, and $\stackrel{\dot{}}{\gamma }$ is the shear rate.

Surface tension is included in this analysis using a potential function as follows:

$P\left(r\right)=\{\begin{array}{ll}{C}_{\text{surface}}{\left(r-r_{\text{e}}\right)}^2\left(\frac{1}{3}r-\frac{2}{3}{r}_{\min }+\frac{1}{2}{r}_{\text{e}}\right)\hfill & r\le r_{\text{e}}\hfill \\ 0\hfill & r>r_{\text{e}}\hfill \end{array}$ (19)

where r_e is the effective radius, C_surface is a constant, and r_min is the initial closest particle distance. As shown in [26], C_surface can be calculated from the surface tension σ. The constant in Equation (19) between the thermoplastic and metallic nozzle or bed is calculated from the contact angle, as shown in [30].

The computational procedures of the MPS method are described in detail in [26]. In this study, not only the pressure but also the non-linear viscosity and heat conduction analysis are implicitly solved using the conjugate gradient method, as shown in a previous study [25]. All material constants are the same as those used in the previous study [25].

3. Experimental Method and 2D Modelling of Nozzle

An FFF-type 3D printer, X-plus, manufactured by Qidi Tech (China), was used for the experiments. Using the 3D printer, a single-path printing was executed first, and the cross-sectional configuration was observed to compare the simulation results. Secondly, a dual-path printing was performed to observe the fusion process between the first and second print paths, as shown in Figure 3. The cross-sectional configurations at the A-a line in Figure 3 were observed, and the configurations of the simulation were compared with the experimental results. To observe the print path configuration, the print path was carefully removed from the bed, and the path was embedded in resin. After polishing the surface of the cross section, the cross-sectional configurations were observed with a video microscope.

A short carbon fibre/PA-6 (Onyx) filament, manufactured by Markforged (USA), was used. The fibre volume fraction was approximately 10%. The nozzle temperature was set to 270˚C. The bed temperature was maintained at 80˚C. Three types of nozzle heights were tested: 0.149, 0.212, and 0.349 mm.

To model a conical 3D nozzle in two dimensions for the MPS model, the plane (y-z plane) orthogonal to the nozzle moving direction (x direction) in Figure 3 was modelled. Figure 4 shows the particle model of a nozzle and bed used for MPS simulations. The cross-section of a nozzle with a hole of 0.4 mm width was used as a 2D nozzle model. As the nozzle tip had a flat area, the flat area and nozzle cylinder were modelled using purple particles, as shown in Figure 4. The other nozzle area, shown as a blue upside-down trapezoid, was omitted in the MPS analysis to reduce the number of particles.

Figure 3. Experimental procedures to confirm the simulation results.

Figure 4. Particle model of a nozzle and bed for MPS analyses.

In almost all FFF-type 3D printers, the filament diameter is 1.75 mm and the nozzle outlet diameter is 0.4 mm. Thus, the specified filament feed speed is not the ejection speed from the nozzle outlet. Moreover, the reaction force of the melted thermoplastic resin from the small-diameter nozzle may cause a slip between the extruder gear and filament [13]. The actual extruded feed speed from the printed path was calculated.

To obtain the exact thermoplastic feed speed v_z from the outlet of the nozzle, pre-experiments were conducted to measure the injected thermoplastic resin cross-sectional areas, and the relationship between the specified filament feed speed v_z and the nozzle movement speed v_x was obtained as follows:

$A_{\text{p}}{v}_{\text{x}}t=\pi {\left(d_{\text{n}}/2\right)}^2{v}_{\text{z}}$ (20)

where A_p is the measured cross-sectional area of the print path, d_n is the nozzle outlet diameter, and t is the printing time. As v_x and t are obtained from the printing condition and A_p can be measured, Equation (20) gives v_z.

As shown in Figure 5 (top view), the nozzle moved in the x-direction during printing. The cross section of the front view was modelled in a 2D analysis, as shown in Figure 5 (front view). The brown particles in Figure 5 (front view) indicate the nozzle. Only the nozzle tip was modelled here. The nozzle center located the 1-1’ cross section at time t₁, as shown in the top view of Figure 5. At time t₂ (t₂ > t₁), the nozzle moved in the positive x-direction by v_xt_d (t_d = t₂ − t₁). As the 2D nozzle was modelled at the center of the nozzle, it indicates the cross-section at the 1-1’ line in the top view of Figure 5.

In the front view of Figure 5, the fused filament diameter d_n is 0.4 mm. The cross-section of the 2D nozzle of the front view was fixed at the coordinate. At time = t₂, the fused filament diameter shrank, as shown by the 2-2’ line in the top view. This was caused by the movement of the nozzle in the x-direction. This shrinkage is shown in the front view of time = t₂ in Figure 5. The outer diameter of the nozzle shrank from 2r_n (r_n is the nozzle wall radius at the nozzle tip) to $2\sqrt{r_{\text{n}}^2-{\left(v_{\text{x}}{t}_{\text{d}}\right)}^2}$, as shown in the top view.

However, it is impossible to remove particles during the MPS analysis. The sudden removal of particles may cause pressure oscillations. The shrinkage of the fused filament in the front view was simulated by a moving nozzle in the y-direction. The movement of the nozzle in the y-direction simulated the shrinking of the fused filament only in the left half of the print. Thus, in this study, the left half of the print was analysed, and the right half of the print was obtained using the relationship of symmetry. This is illustrated in Figure 6.

Let us calculate the nozzle velocity v_y. The distance moved by the nozzle y can be calculated as follows:

$y=r_{\text{n}}-\sqrt{r_{\text{n}}^2-{\left(v_{\text{x}}t\right)}^2}$ (21)

From Equation (21), we can calculate the virtual movement of the nozzle in the y-direction v_y by differentiating Equation (21) with respect to time t:

Figure 5. Three-dimensional nozzle movement on two-dimensional nozzle simulation at the center cross-section. The front view is the two-dimensional model at the cross section on the 1-1’ line, and r_n is the nozzle wall radius at the nozzle tip.

Figure 6. Two-dimensional nozzle movement to simulate three-dimensional nozzle movement.

$v_{\text{y}}=\frac{\partial y}{\partial t}=\frac{v_{\text{x}}^{2}t}{\sqrt{r_{\text{n}}^2-{\left(v_{\text{x}}t\right)}^2}}$ (22)

According to Equation (22), the nozzle moved in the y-direction. However, when v_xt approaches r_p (the center of the nozzle reaches the right edge of the nozzle outlet), the calculated v_y diverges to infinity. To prevent this infinite speed, the nozzle was moved far away before v_xt approached 90% of r_p.

The parameters used in this analysis are shown in Table 1.

The filament used was short carbon/PA-6 Onyx, manufactured by Markforged (USA). The fibre volume fraction of Onyx was 12.1%, and the density was 1180 kg/m³ in the datasheet [31]. The thermal conductivity of Onyx was calculated using the rule of mixtures. The properties of Onyx (short carbon/PA-6) were listed in a previous paper [25].

In the present study, three cases of nozzle height were adopted for the experiments and simulations. The three types of nozzle heights were 0.148 mm (Type A), 0.212 mm (Type B), and 0.329 mm (Type C). The other parameters used in the MPS simulation were the same as those listed in a previous paper [25].

4. Single-Path Experiments and Simulations

Three types of nozzle heights were used to print a single path in the experiments to obtain the cross-sectional area of the print path A_p, which was needed to determine the exact filament seed in the experiments (Equation (20)). The cross sections of the three types of nozzle heights are shown in Figure 7.

From these cross-sectional observations, the cross-sectional area of each print path was measured, and the filament feed speed v_z was obtained using Equation (20). The measured cross-sectional area (A_p) of type A was 0.179 mm². The area of type B was 0.145 mm² and that of type C was 0.151 mm². As the nozzle movement speed v_x was fixed at 18 mm/s in these experiments, the filament feed speed $v_{\text{z}}=4A_{\text{p}}{v}_{\text{x}}/d_{\text{n}}^2=\left(4\times 18/{0.4}^2\right)A_{\text{p}}=450A_{\text{p}}$ mm/s could be calculated. For type A, v_z was 80.55 mm/s. For type B, v_z was 65.25 mm/s, and v_z was 67.95 mm/s for type C.

Table 1. Analytical conditions for calculation by MPS method.

As an example, the MPS simulation results for type B are shown in Figure 8. Figure 8 shows the temperature changes during the printing process. Printing starts at t = 0.0 s, and this is shown in Figure 8(a). Figure 8(b) shows the temperature distribution during printing (t = 0.012 s). Figure 8(c) shows the end of printing (t = 0.0175 s). Figure 8(d) shows the movement of the nozzle (t = 0.0419 s). As mentioned before, the movement of the nozzle on the right side represents the movement of the nozzle in the x-direction.

(a) (b) (c)

Figure 7. Experimental results of single-path printing tests. (a) Cross-sectional view of type A (nozzle height = 0.148 mm); (b) Cross-sectional view of type B (nozzle height = 0.212 mm) (cited from Figure 8 of the reference [25] with permission); (c) Cross-sectional view of type C (nozzle height = 0.329 mm).

(a) (b) (c) (d)

Figure 8. Simulated results of temperature distribution for Onyx obtained by the proposed method based on the MPS method. (a) Distribution of temperature (t = 0.0 s); (b) Distribution of temperature (t = 12 μs); (c) Distribution of temperature (t = 175 μs); (d) Distribution of temperature (t = 419 μs).

Figure 9 shows the results for the liquid fraction. The red particles indicate that the PA-6 resin melted completely. The blue particles indicate that the PA-6 resin was perfectly solid. Figure 9(b) shows that the PA-6 in the particles reaching the bed became solid. This was because the bed temperature was maintained at 80˚C. Figure 9(d) shows that the upper half of the printed composite part was liquid, even after printing was finished.

After solidification, the configuration of the printed path in the 2D simulations was compared with the experiments as shown in Figure 7. In the simulations, the right half was not appropriate because of the nozzle movement. Thus, the right half of the printed part was used, and the entire cross-section was made using symmetry. Smooth circumscribed curves were adopted for the simulation to make the configuration of the printed path cross-section. Similar to the previous research [25], the zero-mean normalized cross correlation (ZNCC) method was used to compare the configurations. The ZNCC method uses a comparison between pixels in grey scale images as follows:

$S_{\text{ZNCC}}\left(x,y\right)=\frac{{\displaystyle \sum {\displaystyle \sum \left\{\left(g\left(x+i,y+j\right)-\bar{g}\right)\left(f\left(i,j\right)-\bar{f}\right)\right\}}}}{\sqrt{{\displaystyle \sum {\displaystyle \sum {\left(g\left(x+i,y+j\right)-\bar{g}\right)}^2}}}\sqrt{{\displaystyle \sum {\displaystyle \sum {\left(f\left(i,j\right)-\bar{f}\right)}^2}}}}$ (23)

(a) (b) (c) (d)

Figure 9. Simulated results of liquid fraction distribution for Onyx obtained by the proposed method based on the MPS method. (a) Distribution of liquid fraction (t = 0.0 s); (b) Distribution of liquid fraction (t = 12 μs); (d) Distribution of liquid fraction (t = 175 μs); (e) Distribution of liquid fraction (t = 175 μs).

In Equation (23), g and f represent the two grey scale values that we want to compare. The upper bar represents the mean value of the images. The calculated SZNCC values are presented in Table 2.

Table 2 shows that the SZNCC values are almost equal to 1. This means that the 2D simulations agreed very well with the experimental results. Therefore, the 2D simulation method was proven to be effective for the 3D printing of short carbon fibre/PA-6 composites.

5. Dual-Path Experiments and Simulations

The same 3D printer was used for the dual-path experiments. Two types of experiments were performed. The four experimental conditions are presented in Table 3. The nozzle height and cross-sectional area were measured from the cross-sectional observations. The path spacing was given the G-code because it was difficult to measure the path spacing from the cross-sectional observations after the fusion process. Figure 10 shows the results of the cross-sectional observations of the dual-path experiments. Figure 11 shows the magnified images of the fused area between the two paths.

The magnified images at the fused area in Figure 11 show that the two paths are completely joined at the path spacing 0.64 mm; however, there are slight gaps at the path spacing 0.70 mm.

Table 2. Correlation coefficient SZNCC between analysis results and experimental results.

Table 3. Two types of dual-path experiments.

(a) (b)

Figure 10. Results of dual-path cross-sectional observations in test specimens. (a) Type D; (b) Type E.

(a) (b)

Figure 11. Magnified images of the fused area between the dual paths. Figures in the parenthesis are given G-code distances (mm) between the dual paths. (a) Type D (0.64); (b) Type E (0.70).

Examples of the fusion simulations are shown in Figure 12. Figure 12 shows the temperature changes during the fusion process for type D. In these simulations, only half of the first path was placed on the left of the nozzle. At t = 419 s, the nozzle moved to the right to represent the shrinking nozzle diameter, as shown in Figure 9.

Figure 12(b) shows that heat conduction occurred from the ejected path to the first path at the interface between the two paths. Figure 12(c) and Figure 12(d) show that the fusion process occurred at the interface. However, there was no high-temperature area in the first path, except for the interface. The resin flowed only around the interface.

To observe the fused configuration between the two paths, smooth circumscribed curves were adopted for this simulation. The results are shown in Figure 13.

When the spacing between the two paths was small (0.64 mm), an almost perfect fusion at the interface was observed in Figure 13(a) (type D). In the experimental result of type D, however, the second path covered the right part of the first path, and the step was made at the top of the fused area. This step was not shown in the simulation results. The difference was caused by the experimental error in the levelling of the print bed. As the print bed was not perfectly horizontal, there was a gap between the nozzle of the second path and the top of the first path. This caused a step in the experimental results. As there was no experimental error in the simulation, there was no step at the top of the fused area.

Figure 13(b) shows the results for type E with a larger spacing (0.70 mm) between the first and second paths. The simulation result showed that there was a small gap at the bottom of the fusion area, and there was no gap at the roof of the fusion area. However, there was a gap at the roof of the fusion area and at the bottom of the fusion area in the experimental result. As the roof of the first and second paths seemed almost at the same height, this did not come from bed levelling. In order to investigate the effect of error of spacing between the paths, we simulated two more cases: the spacings were 0.72 mm and 0.74 mm. The results are presented in Figure 14.

When the spacing was 0.74 mm (Figure 14(b)), the two paths were completely separated. When the spacing was 0.72 mm (Figure 14(a)), the configuration of the fusion area was similar to the experimental result (Figure 13(b)). As the published location error of the printer is 0.011 mm, the difference may be due to the print path location error of the 3D printer. Considering the error of the experiments, the fusion area configuration showed good agreement with the experimental results.

Although the modified MPS method gives excellent results, the limitation of the method should be investigated about the particle radius dimension. As the experimental data of viscosity change caused by fiber volume fraction was limited to 20% [32], further experiments will be required in future.

(a) (b) (c) (d)

Figure 12. Simulated results of temperature distribution of type D using MPS method. (a)t = 0.0 s; (b) t = 140 μs; (c) t = 175 μs; (d) t = 419 μs.

(a) (b)

Figure 13. Comparison of images at the fusion area between the simulation results and the experimental results. The left part is the first path, and the right part is the second path. (a) Cross-sectional comparison of type D (0.64 mm); (b) Cross-sectional comparison of type E (0.70 mm).

(a)(b)

Figure 14. Magnified simulated results of the fusion area between the first path and the second path. The left has a spacing of 0.72 mm. The right has a spacing of 0.74 mm. (a) 0.72 mm; (b) 0.74 mm.

6. Conclusions

Using composite particles made from short carbon fibre/PA-6 composites, the MPS method was adopted to simulate the print path fusion process in an FFF-type 3D printer. To prevent large-scale simulations, a new two-dimensional nozzle movement simulation that enabled the simulation of nozzle movement in an out-of-plane direction was developed. The new method was applied to single-path and dual-path simulations. Experiments were performed to compare the simulation results. The results obtained are as follows:

1) The 3D print two-dimensional simulation method was in excellent agreement with the cross-sectional configuration of the single-path printing of short carbon fibre/PA-6 composites.

2) The method was in good agreement with the experimental results of fusion of two print paths. This method can be used to simulate print path fusion.

Funding

This work was partially supported by a Grant-in-Aid for Scientific Research (c) [grant number 20K05113] from the Japan Society for the Promotion of Science.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Love, L.J., Kunc, V., Rios, O., Duty, C.E., Elliott, A.M., Post, B.K., Smith, R.J. and Blue, C.A. (2014) The Importance of Carbon Fiber to Polymer Additive Manufacturing. Journal of Materials Research, 29, 1893-1898. https://doi.org/10.1557/jmr.2014.212
[2]	Tekinalp, H.L., Kunc, V., Velez-Garcia, G.M., Duty, C.E., Love, L.J., Naskar, A.K., Blue, C.A. and Ozcan, S. (2014) Highly Oriented Carbon Fiber-Polymer Composites via Additive Manufacturing. Composites Science and Technology, 105, 144-150. https://doi.org/10.1016/j.compscitech.2014.10.009
[3]	Postiglione, G., Natale, G.G., Griffini, Levi, M. and Turri, S. (2015) Conductive 3D Microstructures by Direct 3D Printing of Polymer/Carbon Nanotube Nanocomposites via Liquid Deposition Modelling. Composites Part A, 76, 110-114. https://doi.org/10.1016/j.compositesa.2015.05.014
[4]	Matsuzaki, R., Ueda, M., Namiki, M., Jeong, T.-K., Asahara, H., Horiguchi, K., Nakamura, T., Todoroki, A. and Hirano, Y. (2016) Three-Dimensional Printing of Continuous-Fiber Composites by In-Nozzle Impregnation. Scientific Reports, 6, Article No. 23058. https://doi.org/10.1038/srep23058
[5]	Wang, X., Jiang, M., Zhou, Z., Gou, J. and Hui, D. (2017) 3D Printing of Polymer Matrix Composites: A Review and Prospective. Composites Part B, 110, 442-458. https://doi.org/10.1016/j.compositesb.2016.11.034
[6]	Parandoush, P. and Lin, D. (2017) A Review on Additive Manufacturing of Polymer-Fiber Composites. Composite Structures, 182, 36-53. https://doi.org/10.1016/j.compstruct.2017.08.088
[7]	Ning, F., Cong, W., Qiu, J., Wei, J. and Wang S. (2015) Additive Manufacturing of Carbon Fiber Reinforced Thermoplastic Composites Using Fused Deposition Modelling. Composites Part B, 80, 369-378. https://doi.org/10.1016/j.compositesb.2015.06.013
[8]	Ferreira, R.T., Amatte, I.C., Dutra, T.A. and Bürger, D. (2017) Experimental Characterization and Micrography of 3D Printed PLA and PLA Reinforced with Short Carbon Fibers. Composites Part B, 124, 88-100. https://doi.org/10.1016/j.compositesb.2017.05.013
[9]	Ning, F., Cong, W., Hu, Y. and Wang, H. (2017) Additive Manufacturing of Carbon Fiber-Reinforced Plastic Composites Using Fused Deposition Modelling: Effects of Process Parameters on Tensile Properties. Journal of Composite Materials, 51, 451-462. https://doi.org/10.1177/0021998316646169
[10]	Zhang, W., Cotton, C., Sun, J., Heider, D., Gu, B., Sun, B. and Chou, T.W. (2018) Interfacial Bonding Strength of Short Carbon Fiber/Acrylonitrile-Butadiene-Styrene Composites Fabricated by Fused Deposition Modelling. Composites Part B, 137, 51-59. https://doi.org/10.1016/j.compositesb.2017.11.018
[11]	Naranjo-Lozada, J., Ahuett-Garza, H., Orta-Castañón, P., Verbeeten, W.M.H. and Sáiz-González, D. (2019) Tensile Properties and Failure Behavior of Chopped and Continuous Carbon Fiber Composites Produced by Additive Manufacturing. Additive Manufacturing, 26, 227-241. https://doi.org/10.1016/j.addma.2018.12.020
[12]	Yavas, D., Zhang, Z., Liu, Q. and Wu, D. (2021) Fracture Behavior of 3D Printed Carbon Fiber-Reinforced Polymer Composites. Composites Science and Technology, 28, Article ID: 108741. https://doi.org/10.1016/j.compscitech.2021.108741
[13]	Kubota, M., Hayakawa, K. and Todoroki, A. (2022) Effect of Build-Up Orientations and Process Parameters on the Tensile Strength of 3D Printed Short Carbon Fiber/PA-6 Composites. Advanced Composite Materials, 31, 119-136. https://doi.org/10.1080/09243046.2021.1930497
[14]	Talagani, M.R., DorMohammadi, S., Dutton, R., Godines, C., Baid, H., Abdi, F., Kunc, V., Compton, B., Simunovic, S., Duty, C., Love, L., Post, B. and Blue, C. (2015) Numerical Simulation of Big Area Additive Manufacturing (3D Printing) of a Full Size Car. SAMPE Journal, 51, 27-36.
[15]	Makino, M., Fukuzawa, D., Murashima, T., Kawakami, M. and Furukawa, H. (2017) Analysis of Deposition Modelling by Particle Method Simulation. Microsystem Technologies, 23, 1177-1181. https://doi.org/10.1007/s00542-016-3047-4
[16]	Yang, D., Wu, K., Wan, L. and Sheng, Y. (2017) A Particle Element Approach for Modelling the 3D Printing Process of Fibre Reinforced Polymer Composites. Journal of Manufacturing and Materials Processing, 1, 10-21. https://doi.org/10.3390/jmmp1010010
[17]	Bertevas, E., Férec, J., Khoo, B.C., Ausias, G. and Phan-Thien, N. (2018) Smoothed Particle Hydrodynamics (SPH) Modelling of Fiber Orientation in a 3D Printing Process. Physics of Fluids, 30, Article ID: 103103. https://doi.org/10.1063/1.5047088
[18]	Brenken, B., Barocio, E., Favaloro, A., Kunc, V. and Pipes, R.B. (2019) Development and Validation of Extrusion Deposition Additive Manufacturing Process Simulations. Additive Manufacturing, 25, 218-226. https://doi.org/10.1016/j.addma.2018.10.041
[19]	Ouyang, Z., Bertevas, E., Parc, L., Khoo, B.C., Phan-Thien, N., Férec, J. and Ausias, G. (2019) Smoothed Particle Hydrodynamics Simulation of Fiber-Filled Composites in a Non-Isothermal Three-Dimensional Printing Process. Physics of Fluids, 31, Article ID: 123102. https://doi.org/10.1063/1.5130711
[20]	Ouyang, Z., Bertevas, E., Wang, D., Khoo, B.C., et al. (2020) Smoothed Particle Hydrodynamics Study of a Non-Isothermal and Thermally Anisotropic Fused Deposition Modelling Process for a Fiber-Filled Composite. Physics of Fluids, 32, Article ID: 053106. https://doi.org/10.1063/5.0004527
[21]	Zhang, H., Zhang, L., Zhang, H., Wu, J., An, X. and Yang, D. (2021) Fibre Bridging and Nozzle Clogging in 3D Printing of Discontinuous Carbon Fibre-Reinforced Polymer Composites: Coupled CFD-DEM Modelling. The International Journal of Advanced Manufacturing Technology, 117, 3549-3562. https://doi.org/10.1007/s00170-021-07913-7
[22]	Lee, S.W., Cho, S.J. and Kim, W. (2022) Development of Thermo-Fluid Simulation Technique for Extruder and Chamber of FDM-Type 3D Printer for Printing High-Melting-Point Materials. Microsystem Technologies. https://doi.org/10.1007/s00542-022-05270-3
[23]	Farazin, A. and Mohammadimehr, M. (2022) Effect of Different Parameters on the Tensile Properties of Printed Polylactic Acid Samples by FDM: Experimental Design Tested with MDs Simulation. The International Journal of Advanced Manufacturing Technology, 118, 103-118. https://doi.org/10.1007/s00170-021-07330-w
[24]	Fraga Fiho, C.A.D. (2019) Smoothed Particle Hydrodynamics Fundamentals and Basic Applications in Continuum Mechanics. Springer, Cham.
[25]	Imaeda, Y., Todoroki, A., Matsuzaki, R., Ueda, M. and Hirano, Y. (2021) Modified Moving Particle Semi-Implicit Method for 3D Print Process Simulations of Short Carbon Fiber/Polyamide-6 Composites. Composites Part C, 6, Article ID: 100195. https://doi.org/10.1016/j.jcomc.2021.100195
[26]	Koshizuka, S., Shibata, K., Kondo, M. and Matsunaga, T. (2018) Moving Particle Semi-Implicit Method: A Mesh Free Particle Method for Fluid Dynamics. Academic Press, London.
[27]	Voller, V.R. and Prakash, C. (1987) A Fixed-Grid Numerical Modelling Methodology for Convection-Diffusion Mushy Region Phase-Change Problems. International Journal of Heat and Mass Transfer, 30, 1709-1720. https://doi.org/10.1016/0017-9310(87)90317-6
[28]	Tanaka, M. and Matsunaga, T. (2010) Stabilization and Smoothing of Pressure on MPS Method by Quasi-Compressibility. Journal of Computational Physics, 229, 4279-4290. https://doi.org/10.1016/j.jcp.2010.02.011
[29]	Carreau, P.J. (1972) Rheological Equations from Molecular Network Theories. Transactions of the Society of Rheology, 16, 99-127. https://doi.org/10.1122/1.549276
[30]	Kondo, M., Koshizuka, S. and Takimoto, M. (2007) Surface Tension Model Using Interparticle Potential Force in Moving Particle Semi-Implicit Method. Transactions of Japan Society of Computational Engineering and Science, 2007, Article ID: 20070021. (In Japanese) https://doi.org/10.1115/FEDSM2007-37215
[31]	Micro Carbon Fiber Filled Nylon That Forms the Foundation of Markforged Composite Parts. https://markforged.com/materials/plastics/onyx
[32]	Araki, K. and Nakamura, K. (1988) Melt Flow Property of Fiber Reinforced Nylon and Polycarbonate. Journal of the Textile Machinery Society of Japan, 42, 39-47. (In Japanese) https://doi.org/10.4188/transjtmsj.41.3_T39