Flexural Fatigue in a Polymer Matrix Composite Material Reinforced with Continuous Kevlar Fibers Fabricated by Additive Manufacturing

Fatigue bending tests, under controlled displacement, were performed on a polymer matrix composite material reinforced with continuous Kevlar fibers. The samples were fabricated using the Fused Filament Fabrication (FFF) technique in a Markforged Two® 3D printer. The static characterization delivered a flexural modulus of elasticity of 4.73 GPa and flexural strength of 110 MPa. The applied loading corresponded to 92.3, 88.5, 86.2, and 84.7% of the static flexural displacement, giving 15, 248, 460, and 711 cycles for failure. Additionally, two numerical models were created: one using orthotropic properties for static loading conditions; and a second one using isotropic in-bulk properties for fatigue modeling. The second model was able to reproduce the experimental fatigue results. Finally, morphological analysis of the fractured surface revealed fiber breakage, fiber tearing, fiber buckling, matrix cracking, and matrix porosity.


Introduction
Additive Manufacturing (AM) enables the production of geometrically complex parts without post-processing. AM is used in a wide range of materials and allows the reconstruction of layer-by-layer 3D topologies that, with traditional manufacturing processes such as milling, cutting, casting, forging, and welding, would be complicated to manufacture [1][2][3][4][5]. AM could effectively link topological optimization and the final product by delivering final geometries with no additional processing [2]. Furthermore, it has been estimated that printing costs are independent of the batch size [6]. Besides, being such a recent technology, AM is a technique in constant development [2], which has been integrated into manufacturing chains [7]. FFF is by far the most widely used method of AM techniques [2,4,5,7]. The AM technology provided by Markforged ® (Markforged, Watertown, MA, USA) [8] can provide load-bearing parts matching the strength of aluminum [9,10].
On the other hand, fatigue is a common failure mechanism responsible for 90% of rotating equipment failures in the industry [11]. Methods for characterizing materials under fluctuating loads include stress versus life (σ-N) for loads below the elastic limit, strain versus life (E -N) for loads above the elastic limit but below the maximum stress that produce crack nucleation, and the da/dN versus ∆K method when measurable cracks are present, using Paris-type rules to estimate crack growth [11]. In composites, the interlaminar fracture toughness is usually performed [12]. Saleh et al. [13] stated that a composite material is considered to fail when its residual strength is less than 85% of its ultimate tensile

Background
This section gives the baseline knowledge needed to understand the paper.

Additive Manufacturing
AM is a technology created by Hull in 1984 that began with stereolithography, followed by FFF in 1989 [5]. After the Stratasys patent expired in 2010, there was a spike in FFF development. In 2015, Markforged ® (Watertown, MA, USA) [8] obtained a patent to manufacture FFF-printed composite materials reinforced with continuous fibers, giving good dimensional accuracy and stability [7].
FFF is a manufacturing technique in which topology is reconstructed layer by layer from a G-code file that has previously been converted to STL format. [4,5]. This process is done by software in the cloud called Eiger ® for Markforged ® composites [8]. The 1.75 mm polymeric filament, Onyx ® in this case, is fed into an extruder, heated above its melting point (260 • C, approx.), and deposited on a platform where it cools to form the wanted part. The extruder movement is controlled by software that optimizes the movements and alternates matrix and continuous reinforcement deposition. A schematic of the FFF printing process is shown in Figure 1.
Before printing a part, several process parameters must be selected in the Eiger ® (Watertown, MA, USA) web-based software [2]. These are fiber type (carbon, Kevlar, fiberglass, and high-temperature fiberglass), fiber volume fraction, fiber layout type (concentric and isotropic, as shown in Figure 2), matrix fill pattern (rectangular, hexagonal, triangular, as seen in Figure 2, and solid fill), matrix fill density, matrix deposition angle, and fiber deposition angle [15]. In addition, there are new developments in cellular structures that might give improved strength and stiffness over existing ones [31].
In the same way, as in the design with traditional composite materials, the properties of the parts produced by AM depend on the volumetric fiber fraction, the arrangement and angle of the fibers, and the mechanical properties of the reinforcement and matrix [28,29]. Additionally, FFF creates an anisotropy besides what classic composites show [32]. Therefore, it is common to represent its mechanical behavior with an orthotropic model [2] from matrix and fiber properties such as the ones provided by Markforged ® [33], as shown in Table 1. Melenka et al. [34] proposed a model to estimate mechanical properties in AM printed composite materials based on the volumetric average stiffness (VAS) method [35] but accounting for voids in the matrix [36]. Before printing a part, several process parameters must be selected in the Eiger ® (Watertown, MA, USA) web-based software [2]. These are fiber type (carbon, Kevlar, fiberglass, and high-temperature fiberglass), fiber volume fraction, fiber layout type (concentric and isotropic, as shown in Figure 2), matrix fill pattern (rectangular, hexagonal, triangular, as seen in Figure 2, and solid fill), matrix fill density, matrix deposition angle, and fiber deposition angle [15]. In addition, there are new developments in cellular structures that might give improved strength and stiffness over existing ones [31]. In the same way, as in the design with traditional composite materials, the properties of the parts produced by AM depend on the volumetric fiber fraction, the arrangement and angle of the fibers, and the mechanical properties of the reinforcement and matrix [28,29]. Additionally, FFF creates an anisotropy besides what classic composites show [32]. Therefore, it is common to represent its mechanical behavior with an orthotropic model [2] from matrix and fiber properties such as the ones provided by Markforged ® [33], as shown in Table 1. Melenka et al. [34] proposed a model to estimate mechanical  Before printing a part, several process parameters must be selected in the Eiger ® (Watertown, MA, USA) web-based software [2]. These are fiber type (carbon, Kevlar, fiberglass, and high-temperature fiberglass), fiber volume fraction, fiber layout type (concentric and isotropic, as shown in Figure 2), matrix fill pattern (rectangular, hexagonal, triangular, as seen in Figure 2, and solid fill), matrix fill density, matrix deposition angle, and fiber deposition angle [15]. In addition, there are new developments in cellular structures that might give improved strength and stiffness over existing ones [31]. In the same way, as in the design with traditional composite materials, the properties of the parts produced by AM depend on the volumetric fiber fraction, the arrangement and angle of the fibers, and the mechanical properties of the reinforcement and matrix [28,29]. Additionally, FFF creates an anisotropy besides what classic composites show [32]. Therefore, it is common to represent its mechanical behavior with an orthotropic model [2] from matrix and fiber properties such as the ones provided by Markforged ® [33], as shown in Table 1. Melenka et al. [34] proposed a model to estimate mechanical   Table 1, [19] reported 7.8 MPa for average ultimate tensile strength for a triangular filling at 20%. An analytical summary of reported mechanical properties for a wide combination of parameters is available at [9]. Further details on the printing process may be found in [8].

Fatigue
It is commonly accepted that fatigue failure begins with the nucleation of cracks [11]. For the most part, cracks in composites form in layers perpendicular to the load direction [9,26]; this process is called transverse matrix cracking and involves numerous microcracks. Thus, cracks in the matrix are generally the first form of damage in composites and function as a barrier delaying further macroscopic damage [30].
When a structural component is subjected to a repetitive load well below the yield stress, the recommended design method is stress life (σ-N). This behavior is described by the phenomenological Basquin rule [11], as shown in Equation (1).
where σ is stress, the strength reduction rate depends on the number of cycles and materialdependent constants. A straight line is obtained when plotting Equation (1) in a bilogartimic graph (N, σ). However, the tests were done at a deflection inversion ratio, Ry = −1, as shown in Equation (2).
where Y is the beam's deflection. Two studies were found regarding fatigue studies about composite materials printed by Nylon matrix FFF [14,15]; both were done under axial load. This study takes the configuration that showed the best results when combined with different printing parameters [15,16], see Table 2. In [15], the experimental data, which was performed under a load inversion ratio R = −1, was fitted into Basquin's rule [11]. Furthermore, these materials have been shown to exhibit a loss of stiffness under cyclic loads [27,30].

Materials and Methods
The specimens were manufactured according to ASTM D6272-17 (Standard Test Method for Flexural Properties of Unreinforced and Reinforced Plastics and Electrical Materials by Four-Point Bending) for static ASTM D7774-17 (Standard Test Method for Flexural Fatigue Properties of Plastics) for four-point bending fatigue tests. All were printed on a Markforged 2.0 printer using parameters defined in [15,16] with Onyx as a matrix with a 19% of continuous Kevlar fiber as reinforcement, both provided by Markforged. Table 2 shows the parameters used. Finally, an MTS Bionix ® (Minneapolis, MN, USA) servohydraulic 370.02 universal testing machine (UTS) was used for the static and fatigue tests at room temperature. The machine was equipped with a 25 kN load cell.

Static Test
Tension tests were performed using an MTS 634.12f 25 mm (Minneapolis, MN, USA) axial extensometer, as shown in Figure 3a, at a crosshead speed of 0.5 mm/min. Specimen dimensions were chosen according to ASTM D6272-17: 127 mm long, 12.7 mm wide, and 3.2 mm thick. For these tests, end-tabs [37] were used to guarantee the clamp-specimen grip, to reduce the compression stress due to an excessive grip force and, therefore, an unwanted failure due to stress concentration. The static bending tests were made at 6 mm/min crosshead speed, according to ASTM D6272. A dial gauge indicator, shown in Figure 3b, Tension tests were performed using an MTS 634.12f 25 mm (Minneapolis, MN, USA) axial extensometer, as shown in Figure 3a, at a crosshead speed of 0.5 mm/min. Specimen dimensions were chosen according to ASTM D6272-17: 127 mm long, 12.7 mm wide, and 3.2 mm thick. For these tests, end-tabs [37] were used to guarantee the clamp-specimen grip, to reduce the compression stress due to an excessive grip force and, therefore, an unwanted failure due to stress concentration. The static bending tests were made at 6 mm/min crosshead speed, according to ASTM D6272. A dial gauge indicator, shown in Figure 3b, was used in the static bending test to verify the sample's deflection. Finally, the in-house built bidirectional supports are shown in Figure 3c.  where σ is tension, P is force, x is the distance between the lower support and the point of application of the load, b is the width of the specimen, and d is the thickness of the specimen. The ratio for strain, ε, also recommended in ASTM D6272-17, is shown in Equation (4). The experimental stress was obtained according to a relationship recommended in ASTM D6272-17 and shown in Equation (3).
where σ is tension, P is force, x is the distance between the lower support and the point of application of the load, b is the width of the specimen, and d is the thickness of the specimen. The ratio for strain, ε, also recommended in ASTM D6272-17, is shown in Equation (4).
where Y max is the maximum deflection of the beam.

Fatigue Tests
The alternative four-point fatigue tests were done at ambient temperature, under a displacement inversion factor Ry = −1 under controlled displacement applying a 5 Hz sinusoidal load. An example of the control and feedback displacement signals is seen in Figure 4, taken from the MTS suite ® (Minneapolis, MN, USA) control software. For these tests, in-house supports were built to support loads in opposite directions with a 12 mm rod radius of, as shown in Figure 3c. The applied load levels corresponded to 92.3, 88, 86.2, and 84.7% of the maximum deflection obtained from the static tests. Furthermore, Castro and Meggiolaro [11] explained the effect on loading frequency; at higher frequencies fatigue life is shortened by the temperature rise that the polymer matrix is not capable of properly dissipating, most likely due to local cyclic heating. However, that effect is only significant when operating frequencies experiment a magnitude change of 10 or more [38]. Therefore, the results presented here must only be applied within the tested range. rod radius of, as shown in Figure3c. The applied load levels corresponded to 92.3, 88, 86.2, and 84.7% of the maximum deflection obtained from the static tests. Furthermore, Castro and Meggiolaro [11] explained the effect on loading frequency; at higher frequencies fatigue life is shortened by the temperature rise that the polymer matrix is not capable of properly dissipating, most likely due to local cyclic heating. However, that effect is only significant when operating frequencies experiment a magnitude change of 10 or more [38]. Therefore, the results presented here must only be applied within the tested range.

Numerical Modeling
The elastic problem for a composite material can be described using a linear elastic model and the Galerkin formulation to solve the displacements u through the domain Ω, as shown in Equation (5).

Numerical Modeling
The elastic problem for a composite material can be described using a linear elastic model and the Galerkin formulation to solve the displacements u through the domain Ω, as shown in Equation (5).
where ε and E are strain and elastic tensors, respectively, b the body forces, and t the traction vector.
A finite element simulation was carried out on Ansys ® (ANSYS, Canonsburg, PA, USA). Two numerical models were created. The first model used detailed geometry using the ACP (Ansys Composite Module) with orthotropic properties to simulate the static bending test. Figure 5a shows the cross-section of the detailed model as arranged in the ACP model. Figure 5b shows the organization of the different layers of the specimen, and Figure 5c shows the orientation of the triangular fill layers. The second simulation was a simplified model using bulk properties, as suggested in [9]. The simplified model was used to simulate the fatigue tests because fatigue properties were unavailable for each material region.
The properties of the cross-sectional regions for the detailed model were estimated through the Ansys Material Designer plug-in for Ansys Workbench, which is based on a homogenization method through an RVE (Representative Elementary Volume) model. These properties were then verified with the ROM and the Halpin-Tsai equations [39]. The triangular filling region in the simulation was considered a solid material but with equivalent properties calculated based on the equations of cellular solids, see Section 4.2 [35]. For the Onyx ® solid regions (layers and wall), the mechanical properties assigned to the numerical model were obtained based on the model described in [36] and shown in Table 3.
the ACP (Ansys Composite Module) with orthotropic properties to simulate the static bending test. Figure 5a shows the cross-section of the detailed model as arranged in the ACP model. Figure 5b shows the organization of the different layers of the specimen, and Figure 5c shows the orientation of the triangular fill layers. The second simulation was a simplified model using bulk properties, as suggested in [9]. The simplified model was used to simulate the fatigue tests because fatigue properties were unavailable for each material region.   For the fatigue simulation in Ansys ® using the simplified model, the Onyx reinforced with continuous Kevlar fiber specimen was considered a solid beam with orthotropic properties. These properties were estimated analytically through the volume average stiffness method (VAS) proposed in [34]. Table 4 shows the values used.

Microscopy
Micrographs were taken with a scanning electron microscope (SEM) to observe the failed surface. First, samples were plated with gold and then placed into a Vega 3 Tescan SEM working between 5 and 10 keV equipped with a tungsten filament. The gold plating was done in some samples, finding no difference under the SEM with the no plated samples. The need for no plating was attributed to the chopped carbon fiber in the PA6 matrix that made the sample conductive.

Mechanical Properties Estimation
This section details how the mechanical properties used in the numerical models were estimated.

Triangular Fill Properties Estimation
Gibson and Ashby [35] proposed the concept of relative density, p r shown in Equation (6), based on fill density, p x , and fill material density p s , for a triangular fill.
For equilateral triangle fills p r is given by Equation (7), where t is the width of the triangular cell, and l is the length of the triangular cell.
The shear modulus in the different planes (G 12 , G 23 , G 13 ) were calculated with the equations shown in Equation (8) based on Gs, Onyx shear module.
Finally, the Poisson ratio in the different planes (v 12 , v 23 , v 13 ) was estimated by the Equation (9), where v s is Onyx Poisson ratio.

Solid Onyx Region Properties Estimation
The upper and lower walls and layers were characterized following Rodríguez [31], which considered the level of porosity (p 1 ) of the solid regions printed by FFF. Thus, the material is deposited in the Z direction, as shown in Figure 6.

Solid Onyx Region Properties Estimation
The upper and lower walls and layers were characterized following Rodríguez [31], which considered the level of porosity (p1) of the solid regions printed by FFF. Thus, the material is deposited in the Z direction, as shown in Figure 6. The elastic modulus in the longitudinal direction (E1) and in the transversal directions (E2, E3) are given by Equation (10), which is modified by the porosity level p1. According to Papon and Haque [40] p1 is estimated at 10%. The elastic modulus in the longitudinal direction (E 1 ) and in the transversal directions (E 2 , E 3 ) are given by Equation (10), which is modified by the porosity level p 1 . According to Papon and Haque [40] p 1 is estimated at 10%.
For calculating the shear modulus in the different planes (G 12 , G 23 , G 13 ), Equation (11) was used.

Estimation of the Properties for the Continuous Fiber-Reinforced Regions
The continuous Kevlar fiber-reinforced regions, see Figure 5a, were considered a unidirectional ply composed of a Nylon matrix and continuous fiber, with equivalent properties based on the rule of mixtures (ROM) and the Halpin-Tsai equations [39].
The modulus of elasticity in the longitudinal direction (E 1 ) and the transverse direction (E 2 , E 3 ) were calculated with Equation (13). where Vf is fiber content, and ξ is 2 for circular fibers. The shear modulus in the different planes (G 12 , G 23 , G 13 ) were calculated with Equation (14).
where Vf is 0.36 and v m is 0.39. To calculate Poisson's ratios (v 12 , v 23 , v 13 ) Equation (15) was used. Figure 7 shows the results of the uniaxial tension test for the non-reinforced Onyx ® sample performed with an MTS 634.12F extensometer. The sample had the same printing parameters as the bending specimens, as the best performing ones reported in [15], but did not have continuous Kevlar fiber reinforcement. Stress was calculated as average stress; this is force over the initial area, while strain was estimated as elongation over initial length. Figure 7 depicts exemplary results for the tension test, where a 12.5 MPa maximum stress and 10.84% strain at rupture are seen. As a comparison, Barnik et al. [19] reported for solid Onyx a maximum stress of 31 MPa, whereas Markforged listed 40 MPa, see Table 1. For other AM polymeric materials, Parrado [6] reported a maximum of 48.4 ±1.9 MPa for solid PLA but for woven-layered (0 ± 90 • ) samples. Another source of difference could be attributed to the sample s dimensions [9].  Figure 7 shows the results of the uniaxial tension test for the non-reinforced Onyx ® sample performed with an MTS 634.12F extensometer. The sample had the same printing parameters as the bending specimens, as the best performing ones reported in [15], but did not have continuous Kevlar fiber reinforcement. Stress was calculated as average stress; this is force over the initial area, while strain was estimated as elongation over initial length. Figure 7 depicts exemplary results for the tension test, where a 12.5 MPa maximum stress and 10.84% strain at rupture are seen. As a comparison, Barnik et al. [19] reported for solid Onyx a maximum stress of 31 MPa, whereas Markforged listed 40 MPa, see Table 1. For other AM polymeric materials, Parrado [6] reported a maximum of 48.4 ±1.9 MPa for solid PLA but for woven-layered (0 ± 90°) samples. Another source of difference could be attributed to the sample´s dimensions [9]. Moreover, an expression based on the Euler-Bernoulli beam theory showed that the maximum deflection is 1.15 times the deflection measured at the central supports, where the load is applied. The measured deflection was verified with a dial gauge. Figure 8  Moreover, an expression based on the Euler-Bernoulli beam theory showed that the maximum deflection is 1.15 times the deflection measured at the central supports, where the load is applied. The measured deflection was verified with a dial gauge. Figure 8 shows exemplary results of four points bending tests for three samples. An average modulus of elasticity of 4.37 ± 0.65 GPa and average ultimate stress of 109.9 ± 14.4 MPa were obtained.

Fatigue Tests
The tests were carried out under controlled displacement. The specimens were distributed in four displacement levels corresponding to 92.3%, 88.5%, and 86.2% of the maximum flexural displacement, corresponding to ±12, ±11.5, ±11.2, and ±11 mm, respectively. From these tests, the applied stress and number of failure cycles, Nf, were obtained. Results of the tested specimens are shown in Table 5. A drop in stiffness was observed in the load-displacement cycles, as reported by Saleh et al. [13]. It is understood with this behavior that stiffness changes because of the matrix collapse or the rupture of some of the reinforcing fibers. The energy absorbed by the composite material, represented by the shape of the hysteresis loop, seen in Figure 9, results in matrix-fiber separation, generating voids that act as stress concentrators, thus reducing fatigue resistance. This behavior is consistent with the literature where this failure mechanism has been mentioned [26,39]. Exemplary results of this phenomenon are shown in Figure 9, where one can observe that the load required to produce the ± 11 mm constant deflection lowers every loading cycle. It is observed how the maximum force lowers from 170 N at the first cycle, goes to 146 N at ten cycles, 124 N at 20 cycles, and stays about 116 N from 100 all the way to 300 cycles. Moreover, the crack was not visible until beyond the 85% loss of strength predicted by Saleh et al. [13]. Therefore, we found that Saleh et al. [13] criteria is conservative as it predicts shorter lives as the sample experiences.

Fatigue Tests
The tests were carried out under controlled displacement. The specimens were distributed in four displacement levels corresponding to 92.3%, 88.5%, and 86.2% of the maximum flexural displacement, corresponding to ±12, ±11.5, ±11.2, and ±11 mm, respectively. From these tests, the applied stress and number of failure cycles, N f , were obtained. Results of the tested specimens are shown in Table 5. A drop in stiffness was observed in the load-displacement cycles, as reported by Saleh et al. [13]. It is understood with this behavior that stiffness changes because of the matrix collapse or the rupture of some of the reinforcing fibers. The energy absorbed by the composite material, represented by the shape of the hysteresis loop, seen in Figure 9, results in matrix-fiber separation, generating voids that act as stress concentrators, thus reducing fatigue resistance. This behavior is consistent with the literature where this failure mechanism has been mentioned [26,39]. Exemplary results of this phenomenon are shown in Figure 9, where one can observe that the load required to produce the ±11 mm constant deflection lowers every loading cycle. It is observed how the maximum force lowers from 170 N at the first cycle, goes to 146 N at ten cycles, 124 N at 20 cycles, and stays about 116 N from 100 all the way to 300 cycles. Moreover, the crack was not visible until beyond the 85% loss of strength predicted by Saleh et al. [13]. Therefore, we found that Saleh et al. [13] criteria is conservative as it predicts shorter lives as the sample experiences. Polymers 2022, 14, x FOR PEER REVIEW 13 of 20 Figure 9. Exemplary load-displacement loops for a specimen subjected to ± 11 mm for different number of cycles.
Moreover, the force-displacement loops do not seem to behave linearly as the loading and unloading paths do not follow the same route. From 1 to 2, the load raises as the deflection rises, but somewhere in between, the slope raises; the deflection causes a more rigid sample, so, the sample acts as a bending spring. At 2, the maximum load and displacement are reached. The displacement is inverted, and between 2 and 3, the load drops rapidly but at a lower slope than it did between 1 and 2. This may be due to the stored energy in the sample during the loading stroke, loading path 1 to 2. The load goes negative at 4, and the slope raises again between 4 and 5. It repeats the process seen between 1 and 2 but at negative loads. At 5, the sample reaches a peak for both displacement and load. The unloading occurs right after 5 but again at a different slope than between 4 and 5, perhaps in the same manner that it occurred for a positive load (from 2 to 3); the stored energy helps to bring back the sample to zero force. At 6, and very rapidly, the sample reaches almost zero force but is far from the zero-displacement position. Finally, the displacement reaches 1 to start another cycle.

Numerical Simulation
A mesh size analysis was made to guarantee the convergence of the solution. The displacement is not dependent on mesh size, as seen in Figure 10. After the second mesh size, about 17,000 elements, corresponding to an average of 1 mm elements, the displacement remains about the same. With this mesh size, eight load values were analyzed from 0 to 70 N, as seen in Figure 11, where numerical results are plotted and compared with experimental tests. Moreover, the force-displacement loops do not seem to behave linearly as the loading and unloading paths do not follow the same route. From 1 to 2, the load raises as the deflection rises, but somewhere in between, the slope raises; the deflection causes a more rigid sample, so, the sample acts as a bending spring. At 2, the maximum load and displacement are reached. The displacement is inverted, and between 2 and 3, the load drops rapidly but at a lower slope than it did between 1 and 2. This may be due to the stored energy in the sample during the loading stroke, loading path 1 to 2. The load goes negative at 4, and the slope raises again between 4 and 5. It repeats the process seen between 1 and 2 but at negative loads. At 5, the sample reaches a peak for both displacement and load. The unloading occurs right after 5 but again at a different slope than between 4 and 5, perhaps in the same manner that it occurred for a positive load (from 2 to 3); the stored energy helps to bring back the sample to zero force. At 6, and very rapidly, the sample reaches almost zero force but is far from the zero-displacement position. Finally, the displacement reaches 1 to start another cycle.

Numerical Simulation
A mesh size analysis was made to guarantee the convergence of the solution. The displacement is not dependent on mesh size, as seen in Figure 10. After the second mesh size, about 17,000 elements, corresponding to an average of 1 mm elements, the displacement remains about the same. With this mesh size, eight load values were analyzed from 0 to 70 N, as seen in Figure 11, where numerical results are plotted and compared with experimental tests.
Four load levels corresponding to 92.3, 88.5, 86.2, and 84.7% of the flexural displacement were simulated for the alternative fatigue results. The fatigue life of the material was compared between experimental results and the numerical model. The fatigue simulation used the simplified numerical model described in Section 3.3. The results are seen in Figure 12, where both plots (numerical and experimental) run parallel and very close to each other, with a low as 2.5% difference.   Additionally, results from similar samples (Kevlar reinforced and triangular filling, but axially loaded instead of bending and PA 6 matrix instead of Onyx ® ) retrieved from the literature [15] are presented in Figure 12. In an axially loaded rectangular sample, the normal stress is assumed as average in the whole section, whereas in a bending sample, only half of the section is in tension and the other half in compression. However, in bending, the stress gradient depends on the distance from the neutral plane. The bending samples show lower stress than the axially loaded counterparts. Thus, besides the stress gradient, this stress difference may be attributed to the short fiber embedded in the Onyx ® matrix, present only in bending samples. The space left on the Onyx ® matrix for the chopped carbon fiber may act as a stress concentration factor lowering the fatigue endurance.   Additionally, results from similar samples (Kevlar reinforced and triangular filling, but axially loaded instead of bending and PA 6 matrix instead of Onyx ® ) retrieved from the literature [15] are presented in Figure 12. In an axially loaded rectangular sample, the normal stress is assumed as average in the whole section, whereas in a bending sample, only half of the section is in tension and the other half in compression. However, in bending, the stress gradient depends on the distance from the neutral plane. The bending samples show lower stress than the axially loaded counterparts. Thus, besides the stress gradient, this stress difference may be attributed to the short fiber embedded in the Onyx ® matrix, present only in bending samples. The space left on the Onyx ® matrix for the chopped carbon fiber may act as a stress concentration factor lowering the fatigue endurance. So, experimental data from Figure 12 were adjusted to Basquin's rule [11] Equation (1), giving the results shown in Equation (16) that represent 95.2% of the experimental data. However, the simulation showed a trend to overpredict failure cycles for the upper part of the simulated cycles. For example, the observed cycles for the +/−11 mm were 711, whereas the simulation gave 755. On the other hand, the value of parameter b in Equation (16) was adjusted to −0.074, showing that it is indeed a brittle composite, as previously reported [15]. Additionally, results from similar samples (Kevlar reinforced and triangular filling, but axially loaded instead of bending and PA 6 matrix instead of Onyx ® ) retrieved from the literature [15] are presented in Figure 12. In an axially loaded rectangular sample, the normal stress is assumed as average in the whole section, whereas in a bending sample, only half of the section is in tension and the other half in compression. However, in bending, the stress gradient depends on the distance from the neutral plane. The bending samples show lower stress than the axially loaded counterparts. Thus, besides the stress gradient, this stress difference may be attributed to the short fiber embedded in the Onyx ® matrix, present only in bending samples. The space left on the Onyx ® matrix for the chopped carbon fiber may act as a stress concentration factor lowering the fatigue endurance.
So, experimental data from Figure 12 were adjusted to Basquin's rule [11] Equation (1), giving the results shown in Equation (16) that represent 95.2% of the experimental data. However, the simulation showed a trend to overpredict failure cycles for the upper part of the simulated cycles. For example, the observed cycles for the +/−11 mm were 711, whereas the simulation gave 755. On the other hand, the value of parameter b in Equation (16) was adjusted to −0.074, showing that it is indeed a brittle composite, as previously reported [15].
Finally, the difference between experimental and simulated results went from 33.3% at low cycles to 6.2% for the longest recorded cycles. Figure 13 shows a general view of three fatigue surfaces taken by a digital camera. The macroscopic crack is nucleated in the polymer matrix on the most strained surface and propagates until it reaches the fiber layer. This nucleation can be attributed to the pure bending moment in the sample s central zone, which produces normal stress. The combination of the remaining Onyx ® matrix, fibers, and nucleated surface cracks account for the elasticity module of the composite. The Kevlar fibers, which are stronger but more brittle, fracture first, starting with the closest to the external surface and consequently changing the composite stiffness. Thus, the elastic modulus changes with every fractured fiber. Such behavior explains the observed decrease in absorbed energy seen in Figure 9. Some samples presented wrinkles on the surface with no evident crack, see Figure 13. Such morphology may be attributed to fiber buckling [9,30]. The continuous fibers are slender, and during the reversal load, they experience compression, which may induce buckling, whereas the matrix has not failed yet. Such buckling may induce fiber-matrix separation. Therefore, composites may fail due to fatigue even under a negative load.

Failure Analysis
Samples were submerged in liquid nitrogen and fractured by a sudden load to observe the failed surface. Figure 14 shows layers of molded polymer (Onyx), with their respective aggregates of short carbon fiber distributed in the Nylon matrix; the distributed fibers appear to be stress concentrators. Empty spaces can be seen between each layer of polymer, and details of the porosity can also be observed in the polymeric matrix. The successive passes done to deposit the molted polymer left void spaces, as documented by Rodríguez [36]. Ekoi et al. [17] also reported a high porosity. Therefore, a high fiber content value may reduce the composite's strength because of the higher voids between matrix and fiber [9]. In this study, the fiber content was fixed, but it is a parameter to consider when designing a composite part. Figure 14 also shows the reinforcing fibers on each side of the Onyx ® matrix. An enormous difference in the fracture shape is observed, the reinforcing fiber is thinner, and its fracture shape is slipped. Figure 15a shows a fiber pullout from the matrix in an axial static test; mechanism reported in literature [17,23,30]. While Figure 15b shows fiber tearing because of the alternative bending fatigue test stress, in Figure 15a we can conclude that applied load is transmitted to the matrix, but the fibers are the ones withstanding the load. Such arrangement creates a strain gradient that pulls the fiber out of the matrix when it reaches a limit. Figure 15b shows the fiber for an alternative bending sample. We can observe the failed surfaces revealed fiber breakage mechanism, fiber tearing, and fiber buckling. Such mechanisms were summarized by Awaja [30]. Samples were submerged in liquid nitrogen and fractured by a sudden load to observe the failed surface. Figure 14 shows layers of molded polymer (Onyx), with their respective aggregates of short carbon fiber distributed in the Nylon matrix; the distributed fibers appear to be stress concentrators. Empty spaces can be seen between each layer of polymer, and details of the porosity can also be observed in the polymeric matrix. The successive passes done to deposit the molted polymer left void spaces, as documented by Rodríguez [36]. Ekoi et al. [17] also reported a high porosity. Therefore, a high fiber content value may reduce the composite's strength because of the higher voids between matrix and fiber [9]. In this study, the fiber content was fixed, but it is a parameter to consider when designing a composite part. Figure 14 also shows the reinforcing fibers on each side of the Onyx ® matrix. An enormous difference in the fracture shape is observed, the reinforcing fiber is thinner, and its fracture shape is slipped. Figure 15a shows a fiber pullout from the matrix in an axial static test; mechanism reported in literature [17,23,30]. While Figure 15b shows fiber tearing because of the alternative bending fatigue test stress, in Figure 15a we can conclude that applied load is transmitted to the matrix, but the fibers are the ones withstanding the load. Such arrangement creates a strain gradient that pulls the fiber out of the matrix when it reaches a limit. Figure  15b shows the fiber for an alternative bending sample. We can observe the failed surfaces revealed fiber breakage mechanism, fiber tearing, and fiber buckling. Such mechanisms were summarized by Awaja [30].

Conclusions
The alternative bending fatigue behavior of a polymeric matrix composite material reinforced with continuous Kevlar fiber and printed by the novel Fused Filament Fabrication technique was analyzed. For this, the geometry of the specimens was established

Conclusions
The alternative bending fatigue behavior of a polymeric matrix composite material reinforced with continuous Kevlar fiber and printed by the novel Fused Filament Fabrication technique was analyzed. For this, the geometry of the specimens was established according to the ASTM D6272-17 standard. Static testing of said specimens gave a flexural modulus of elasticity of 4.73 GPa and flexural strength of 110 MPa.
The experimental stress versus the number of cycles (S-N) curve was obtained for the composite material. Alternative bending tests determined the curve according to ASTM D7774-17. From these tests, the failure cycles of the specimens were found, 15, 248, 460, and 711 cycles, corresponding to 92.3, 88.5, 86.2, and 84.7% of the maximum displacement, respectively. Finally, the composite material stiffness degradation was evidenced by the application of cyclic loads.
Numerical models of the manufactured specimens were created, and the static and fatigue behavior of the material was simulated using the ANSYS composite module. Results were obtained for mechanical properties such as modulus of elasticity and flexural strength, 4.48 GPa and 113 MPa, respectively. In the same way, numerically simulated cycles were very close to those obtained experimentally. That is 20, 250, 462, and 755 failure cycles for 92.3, 88.5, 86.2, and 84.7% of the maximum displacement, respectively. Thus, these numerical results were validated against experimental tests. Consequently, the computational model was validated using bulk mechanical properties instead of a more complex model that would require mechanical fatigue properties for each component of the composite material.
Inspection of the failed surfaces revealed the mechanisms of fiber breakage, matrix cracking, and matrix porosity for the static tests, whereas for alternative tests fiber tearing, fiber buckling, matrix cracking, and matrix porosity. The matrix porosity is seen as empty spaces between adjacent polymer layers. These types of failure mechanisms agree with the literature. Matrix cracking started at the site of maximum normal stress. No evidence of interlayer slipping or printing defects, other than commonly reported porosity, were observed.