Mechanism based failure of 3D-printed continuous carbon fiber reinforced thermoplastic composites

The present work describes a computational mechanism based failure analysis conducted for 3D-printed continuous carbon fiber reinforced thermoplastic composites (CFRTPCs), which could not be seen in the available literature. The material failure is investigated based on intraply failure evaluation and adopts different failure criteria for the material constituents. The micromechanical modeling employs the Asymptotic Homoge-nization technique and comprises the selection of a representative volume element statistically equivalent to the microstructure of the material, which is identified from cross-section micrographs. In contrast to recent work, it is demonstrated that an additional relation is required for the macroscopic deviatoric stresses acting over the matrix. This avoids an overestimation of the matrix failure when the reinforced lamina is subjected to longitudinal and shear loads. The resulting failure envelopes are presented and compared to those provided by analytical failure theories available in the literature. The results obtained by the micromechanical approach showed its ability to predict failure of 3D-printed CFRTPCs, in addition to bring different elements for the discussion that could not be captured with analytical models. In this context, it is believed that the characteristics inherent to the microstructure reproduced in the RVE, particularly contributed to obtaining more realistic failure envelopes.


Introduction
Nowadays, the material extrusion based technology Fused Filament Fabrication (FFF) [1][2][3][4][5], which was initially designed to work with unreinforced thermoplastic materials, is currently able to produce fiber reinforced thermoplastic parts.In general, the FFF process is also referred by the term 3D-printing, which has been widely employed in scientific investigations.3D-printed composites reinforced by continuous fibers have potential to become important options in the production chain of structural components, due to increased production flexibility and expected material enhancements in view of current research efforts.According to the available literature, there are several studies investigating their elastic mechanical properties employing both numerical and experimental techniques [6][7][8][9][10][11][12][13][14][15][16][17][18].Over the last years, most of the conducted research has been supported by the Markforged modified FFF process [19], although conventional 3D-printers have been also adapted to work with continuous fiber reinforced filaments [20][21][22][23][24][25].The Markforged printers MarkOne® and MarkTwo® are able to produce parts reinforced by continuous carbon fibers, Kevlar® or fiberglass where the filament is embedded into a thermoplastic matrix prior to nozzle extrusion.Additional information about this modified FFF process are available in [15,19,26].In despite of the recent valuable contributions, there is still a lack of studies concerning the failure of 3D-printed composite materials.
In order to better understand the failure of unidirectional fiber reinforced 3D-printed composites, it is first necessary to investigate the failure mechanism acting over its constituents.In general, traditional fiber reinforced composite materials can undergo several types of failure at the constituent level, and continuous fiber reinforced 3D-printed composites are included in this group.In recent years, several investigations on the mechanism based failure of traditional fiber reinforced composite materials have emerged.Among the mechanism based failures, interfacial debonding, brittle matrix cavitation, and dilatational failure in ductile matrices are of high relevance in the investigations [27][28][29][30][31].It is also known from the literature that the material characteristics at the constituent level, such as the fiber arrangement, consequently affect the stresses at the microscopic level [32][33][34][35][36]. Furthermore, lamina failure may have different mechanisms depending on where the critical points exist.Therefore, appropriate numerical models and failure criteria for each constituent should be used to predict the lamina failure.
According to the literature survey, it is possible to verify the importance of the micromechanical analysis in the prediction of failure in fiber reinforced composite materials [37][38][39][40][41][42][43][44][45][46][47][48].For instance, phenomenological failure theories with physical considerations still have limitations since the laminated composite material is considered homogeneous and orthotropic.In fact, these approaches should consider the material as heterogeneous and anisotropic instead.Besides, assuming a laminated composite material free of voids and unexpected inclusions, typically the failure process initiates on the matrix or interface.Specifically, defects inherent to the manufacturing process may affect the failure of these constituents, which are not considered in the phenomenological approaches.
Among the difficulties in predicting constituent failure throughout micromechanical analyses, a lack of information exists concerning the correlation between stresses at macroscopic level and stresses at microscopic level [39].In this context, the use of Asymptotic Homogenization [49][50][51][52] plays an important role.This technique is able to predict the homogenized equivalent properties of the heterogeneous material and also determine the stresses at the microscopic level.For this purpose, the macro and microstructural responses are obtained from two uncoupled problems, one in each scale.Then, both responses can be superimposed in order to identify the influence of the heterogeneity on the macro-structural response which circumvents the difficulty noted in [39].For instance, the Asymptotic Homogenization technique was successfully applied in the construction of failure envelopes for traditional fiber reinforced composite materials [44].It is worth noting that the mathematically rigorous approach of Asymptotic Homogenization eventually leads to a complex in-house implementation, mainly when complicated microstructures are considered.In general, only classical fiber arrangements, e.g.square and hexagonal, have been considered.In order to overcome these difficulties, a complete methodology to implement the Asymptotic Homogenization using ABAQUS® was presented in [53] and includes the assessment of the stresses at constituent level.

Objectives and contributions
In view of the lack of studies in the literature, the present work is aimed at conducting micromechanical analyses in order to predict the failure of the 3D-printed continuous carbon fiber reinforced thermoplastic composites, hereinafter referred to as 3D-printed CFRTPCs, supplied by Markforged [19].To this end, a procedure is described for analyzing failure at the constituent level adopting different failure criteria for the material constituents (see Section 2) and taking into account a representative volume element (RVE) that has a microstructure statistically equivalent to the microstructure of the material.The statistically equivalent RVE, obtained according to statistical spacial descriptors, is presented in Section 3 and its accuracy is verified comparing the homogenized properties with experimental data available in the literature.The numerical homogenization models are based on the Asymptotic Homogenization technique, whose implementation followed the methodology described in [53].In contrast to recent work, a different approach is presented in Section 4 to determine both the matrix and interface strength properties.The procedure adopted for the construction of the failure envelopes is presented in Section 5 and a special condition for the matrix failure, under transverse tension and compression, is proposed to better predict the failure envelopes of the 3D-printed CFRTPCs.In particular, this special condition for the matrix failure could not be seen in the available literature.Lastly, the micromechanical failure envelopes for 3D-printed CFRTPCs are presented in Section 6 and compared to those provided by analytical failure theories.Finally, the discussion about the obtained results is supported by a geometrical comparison between the failure envelopes.

Fiber failure criterion
A certain stress state is usually regarded to be the limit for fiber failure in a unidirectional fiber reinforced composite lamina when, under combined loading, the macroscopic stresses σ0 11 , along the fiber direction, reach the longitudinal strength, determined by uniaxial testing.From the micromechanics of failure point of view, a failure criteria related to the fiber is then required.Since the fibers are connected along transverse directions only through the matrix, the failure criterion in tension FFC T , and in compression FFC C , can be simplified to the maximum longitudinal stress failure conditions as [39]. (2)

Matrix failure criterion
Since matrices are usually more sensitive to tensile stresses than to compressive stresses, Ha et al. [39] adopted an expression for the matrix failure criterion in terms of the stress tensor invariants at the microscopic level in the matrix, which can be written as where the von Mises equivalent stress σ VM is given by the first invariant of the stress tensor I 1 is the critical value for the deviatoric stress invariant σ cr VM is the critical value for the volumetric stress invariant I cr 1 is α = C m /T m and η m is a power parameter.T m and C m are respectively the tensile and compressive strength of matrix.When η m = 2, Eq.(3) can be considered as a particular case of the tensorial failure criterion proposed by Tsai and Wu [54] when applied for isotropic materials with different tensile and compressive strengths.A closer inspection of Eq. (3) reveals that it might be considered as a yield criterion instead of a failure criterion.Moreover, it is particularly similar to the well known Drucker-Prager yield criterion [55].As a consequence, the results obtained with Eq. ( 3) are considerably conservative, and this is more pronounced when applied to ductile matrices.Nevertheless, the present work adopted this approach, since it has been successfully applied in [39,44].

Interface failure criterion
The interface between fiber and matrix may fail due to normal and tangential tractions.The quadratic failure criterion adopted in [39,44] considers both effects on the interaction between them according to where t n is the interfacial traction along the normal direction, t s is the interfacial shear traction, Y n is the interfacial traction strength along the normal direction, and Y s is the interfacial shear strength.The angular brackets 〈〉 returns the argument if positive and zero otherwise.It can be seen in the available literature that interfacial debonding may be considered as an initial failure mechanism of polymer matrices in fiber reinforced composite materials [30,31].Based on a conservative approach, the present work assigned material failure when the failure condition in Eq. ( 8) is satisfied.A similar approach was adopted in [44].

Numerical modeling
As previously mentioned, the characteristics at the constituent level, such as fiber arrangement, consequently affect the stress distribution at the microscopic level.Different cross-section micrographs of the 3Dprinted CFRTPCs herein investigated are displayed in [15].From those images, it is possible to verify that the fibers have a substantial random placement, instead of being arranged in a regular array, such as square or hexagonal.This characteristic was also identified previously [13,56] and can be associated with a low fiber volume fraction and/or low compaction.In addition, some degree of agglomeration or clustering is also observed.Besides, it can be verified from the cross-section micrographs displayed in [15] that fibers are mostly circular in cross-section.
Taking into account the characteristics at microscopic level mentioned above, a computationally feasible RVE required to model the behavior of the 3D-printed CFRTPCs shall have its constituents in a spatial distribution statistically equivalent to that found in the microstructure of the whole domain.In other words, it means that a region should not be randomly cropped from the cross-section micrographs under the risk of not being representative of the whole microstructure in terms of agglomeration, periodicity and manufacturing-induced characteristics.Therefore, to avoid these issues, the present work applied a methodology to select an equivalent microstructure from the main crosssection micrographs by cropping several regions with different sizes and from different locations, which were characterized according to different spatial descriptor functions.Detailed information about the adopted methodology is given as follows.
From recent literature [35,[57][58][59][60][61][62][63][64] it can be verified that several statistical techniques have been applied for characterizing spatial distributions of individuals in populations.Among them, the calculation of the nearest neighbor distances and the second order intensity function (also known as Ripley's K-Function [65]) have been widely applied on the generation of RVE based models of continuous fiber reinforced composite materials [35,57,58,63].In view of this reason, these techniques were adopted in the present work.In a set of points arranged in a spatial distribution, the nearest neighbor distance can be defined as the smallest Euclidean distance between a reference point and its neighbors.

Thus, defining P
} as a set of n p points containing the centers of fibers, distributed in the cross-sectional area of an unidirectional fiber reinforced composite, the nearest neighbor distance di between a fiber p i ∈ P, with center (x i , y i ), and its neighboring fibers p j ∈ P, with centers (x j , y j ), is In contrast to recent work, the mean and the standard deviation values of the nearest neighbor distances were analyzed separately.Defining r f as the radius of fibers in a distribution, it can be verified that the mean value μ( d /r f ) tends to 2 if the fibers are in a condition of maximum agglomeration, e.g. in a hexagonal arrangement.Higher values of μ( d /r f ) indicate that the fibers in the distribution are more spaced or dispersed.The standard deviation of the nearest neighbor distances σ( d) were computed in order to provide information about periodicity.If σ( d) tends to 0, it can be verified that the fibers are arranged in a periodic distribution, since all the nearest neighbor distances are equal.The more distant from zero, the less periodic is the distribution.
Taking into account the set of n p points P previously defined, in addition to a set of line segments delimiting the region of study L = {l 1 , l 2 , …, l n } and the area of study A, for a given radial distance h k and center points p ∈ P the second order intensity function K(p, h k , L) can be defined as where ) is the Euclidean distance between the centers of fibers p i (x i , y i ) and p j (x j ,y j ), and w(p i , h k , L) is a weight function that takes into account the edge effects returning the proportion of the circumference with radius h k contained within the region of study bounded by the line segments L to the whole circumference with radius h k .Referring to In regards to the characterization of fibers distributed in the crosssectional area of an unidirectional fiber reinforced composites, if the plot K(h) × h provides a stair-shaped response it means that the fibers are distributed in periodic arrangements, e.g.square or hexagonal.Comparing the plot of K(h) to the plot K p (h) it is possible to determine if the fibers are either agglomerated or dispersed.If K(h) > K p (h), then the fibers in the distribution are more agglomerated (or clustered).Contrariwise, if K(h) < K p (h) the distribution presents some degree of dispersion, and consequently, the fibers are more spaced.
After selecting suitable spatial descriptors for the statistical characterization, and defining their respective conditions and criteria, a bank of images was generated based on cropped regions of three cross-section micrographs from different samples of the 3D-printed CFRTPCs.The size of the cropped regions were selected according to a parameter δ, defined as a multiplication factor applied over the fiber radius where the edge sizes of the cropped regions are given by l edge = r fiber × δ.More details can be seen in [35].A reference value for the fiber radius r fiber , was initially measured with the support of scanning electron microscope images.In this case, a fiber radius r fiber ≈ 4μm was found.This value was later confirmed analyzing the cross-section micrographs of the 3D-printed CFRTPCs.In order to cover a wide range of sizes of cropped regions, i.e. from a minimum representative size up to the limit imposed by the vertical size of the cross-section micrograph, the parameter δ was set to δ = {10, 20, 30, …, 150}.For each parameter δ, a set of ten non-intersecting random regions was defined.Therefore, in the present work 150 regions were cropped from each one of the cross-section micrographs totalizing 450 analyzed images.
According to the periodicity criteria adopted for nearest neighbor distance, cropped regions with size δ ≥ 50 are representative of their respective main regions of the cross-section micrograph.Additionally, according to the agglomeration criteria adopted for same technique, it is suggested that cropped regions with size δ = 50 or δ = 60 are suitable to represent their respective main regions of the cross-section micrograph.With respect to the results obtained for the second order intensity function, those obtained for δ ≥ 40 were able to represent the main region of their respective cross-section micrographs in terms of periodicity and agglomeration.Table 1 summarizes the minimum recommended size δ according to the adopted criteria in the statistical characterization.Taking into account the overall results in Table 1, it can be verified that cropped regions with size δ ≥ 50 are statistically representative of their respective main regions of the cross-section micrograph.These results are in agreement to those obtained in [35] where a minimum recommended size of a RVE, representing a carbon fiber reinforced epoxy, is when δ = 50.
In an attempt to verify the accuracy of the predicted homogenized elastic properties, the influence of the representative volume element size, and consequently the statistically equivalent microstructure, was then investigated.The numerical homogenization models were based on the Asymptotic Homogenization technique whose implementation followed the methodology fully described in [53].The models were implemented in ABAQUS® software with the application of periodic boundary conditions, which are based on the assumption that the volume element is part of a much larger domain composed by repetitive volumes of identical shape and content.Table 2 presents the mechanical properties of the carbon fiber and resin matrix [6,7] adopted in the finite element discretization applied to compute the homogenized elastic properties.
The finite element model used eight-node linear hexahedral elements with reduced integration of type C3D8R available on ABAQUS®.More details about the element formulation can be found in [66].As an example of the finite element models adopted in the current analysis, Fig. 1 displays the finite element discretization applied to a cropped region with size δ = 50.In Fig. 1 the cropped region was discretized into 63001 eight-node linear hexahedral elements of type C3D8R.It is worth noting that the fibers were assumed as continuous in the finite element discretization, leading the solution obtained by the Asymptotic Homogenization method to be independent of the coordinate along the fiber direction.Consequently, the mesh was not refined along the fiber direction in the finite element models.
The homogenized elastic properties E H 1 , E H 2 , and G H 12 obtained from the RVEs were compared in detail with those experimentally measured in [15].Fig. 2 shows the comparison.For these analyses, RVEs corresponding to cropped regions with size δ = [10,50] were employed.It is very important to re-emphasize that 30 cropped regions, from three different cross-section micrographs, were obtained for each one of the sizes δ = [10,50].Among them, one cropped region for each size δ was selected according to their results, i.e. those whose results were the closest to the result obtained for the main region of their respective cross-section micrograph, with respect to the adopted statistical criteria.It can be verified from Fig. 2 that the differences between the homogenized properties and experimental data computed for E 1 are negligible for all RVE sizes.The differences computed for E 2 and G 12 are not negligible though they are relatively small for the RVE with size δ = 50, when compared to the results obtained for the other RVE sizes.
Taking into account the results from the statistical characterization, described in Table 1, in addition to those in Fig. 2, it was concluded that the RVE with size δ = 50 was able to reproduce, with strong fidelity, both the fiber spatial distribution and the homogenized elastic properties of the 3D-printed CFRTPCs herein investigated.Both characteristics are of great importance to the mechanism based failure analyses.Therefore, the model previously presented in Fig. 1 was adopted in the computation of the micromechanical failure envelopes.In an attempt to provide more details, the homogenized elastic properties obtained for the different sizes δ, are shown in Fig. 3. Typically, orthotropic laminae of traditional composite materials are also considered as transversely isotropic.However, the same hypothesis cannot be adopted for the 3Dprinted CFRTPCs since the plane 2-3 is not a plane of isotropy, as seen in the cross-section micrographs.From Fig. 3 it is possible to verify that the differences between the values of predicted elastic moduli E 2 and E 3 are substantially small, but not zero.Similarly, the differences between the values of predicted shear elastic moduli G 12 and G 13 are also small but not zero, which corroborates the previous statement.

Determination of constituents strengths in 3D-printed CFRTPCs
Since the strength properties of the individual constituents in the 3Dprinted CFRTPCs were not available, a methodology to numerically predict them was adopted.This methodology is based on the procedure presented in [39], which includes a failure criterion fitting that minimizes the orthogonal distance between a cloud of points and its respective failure curve.The methodology follows the algorithm described in Fig. 4(a).Although the strengths of the matrix and the fiber could be alternatively determined from experimental investigations, it is not straightforward to obtain the interface strength properties from experimental measurements [27][28][29][30] which eventually leads to the application of numerical approaches for this purpose [31,36,40].
In regards to the numerical computation of constituents strengths, firstly macroscopic stresses        numerical model.From the Asymptotic Homogenization, the stresses at microscopic level are computed for all the elements in the discretization.Then, using the constituent failure criteria defined in Section 2, the constants T f , C f , T m , C m , Y n , and Y s are determined.
The determination of the longitudinal tensile and compressive strengths of fibers, based on the defined fiber failure criterion, is simply direct as can be seen in Section 4.1.From this point, the procedure adopted for predicting the interfacial tensile strength along the normal direction and the interfacial shear strength is similarly straightforward.The details are presented in Section 4.3.However, the determination of the matrix tensile and compressive strengths involves a more complex process.In their work, Ha et al. [39] assumed the power parameter η m in Eq. ( 3) to be η m = 2.Although their assumption substantially simplified the determination of the matrix strengths, the resulting failure envelope for the matrix in the glass/epoxy composite material investigated in their work was clearly overestimated with respect to the critical value of the first stress tensor invariant.In order to avoid this overestimation, the present work adopted a strategy slightly different from that, where an optimization technique was employed to find the strength parameters which minimize the matrix failure envelope.In the next sections, the constituent strengths of 3D-printed CFRTPCs are determined and details about the adopted procedure are given as well.

Fiber strengths
The fiber strengths T f and C f were determined by analyzing the resulting stresses σ F 11 along the fiber direction at microscopic level for all of the elements representing the fibers.Therefore, a cloud of points defined by the set with n f the total number of elements representing the fiber, was plotted on the space σ 11 × elements.Thus, taking into account Eqs. ( 1) and ( 2), both the tensile and compressive strengths of the fiber are Fig. 5 shows the cloud of points σ F 11 plotted for each element in the discretized domain.As already expected, the fiber tensile strength T f is obtained when macroscopic stresses { σ0 } equivalent to the uniaxial tensile strength of the lamina X T are applied to the RVE.Analogously, the fiber compressive strength C f is obtained when macroscopic stresses equivalent to the uniaxial compressive strength of the lamina X C are applied.The cloud of points obtained when applying macroscopic stresses { σ0 } equivalent to Y T , Y C and S 12 were considerably lower than those obtained when applying macroscopic stresses equivalent to X T and X C .For this reason, they are not shown in Fig. 5 nor were considered in the computation of fiber strengths.From the results shown in Fig. 5, the tensile and compressive strengths of the continuous fibers reinforcing thermoplastic in the 3D-printed CFRTPCs herein investigated, are found to be respectively T f = 2490 MPa and C f = 1558 MPa.

Matrix strengths
In order to determine the matrix strengths, the macroscopic failure stresses X T , X C , Y T and Y C were applied to the RVE, the deviatoric stress invariant σ VM , and the volumetric stress invariant I 1 were computed for all the elements corresponding to the matrix in the discretized finite element model.In this case, the macroscopic failure stress equivalent to the longitudinal shear strength S 12 was not taken into account in the determination of the matrix strengths.This strategy is based on the assumption adopted in [39], in which the ply failure under longitudinal shear does not occur due to the initial constituent failure but occurs due to the cumulative effect after initial failure, e.g.propagation of micro cracks in the matrix.Then, similar to the fibers, an initial cloud of points is plotted on the space I 1 × σ VM .Thus, let M be a set of n B points located at the boundary of the cloud of points plotted in the space Therefore, taking into account the points ( at the boundary of the cloud of points and the matrix failure criterion defined in Eq. ( 3), the function that returns the expected value for the volumetric stress for a given point i, can be written as Let d be the total difference between the expected value for the volumetric stress invariant and the respective stress invariant located at the boundary of the cloud of points such as Analyzing Eq. ( 17), it is possible to realize that, as the difference decreases, the matrix failure criterion curve is closer to the VM .Therefore, an optimization problem can be proposed in order to find the matrix strength parameters T m , α, and η m that minimizes the difference d In the present work, the Simulated Annealing method was applied in order to solve the optimization problem proposed in Eq. (18).More details about the Simulated Annealing method, including its workflow, can be seen in [67].The matrix strength parameters T m , α, and η m were considered discrete variables and the search intervals were defined according to the information available in the literature.For the matrix tensile strength T m , the search interval was defined according to the range mentioned by Markforged in [19].With respect to the parameter α, which gives the ratio between the compressive and tensile matrix strength, the search interval was adopted according to the typical range for neat resin described in [39].Lastly, the search interval for the parameter η m was defined in order to cover a wide range of influence on Eq. ( 3).Thus, the search intervals were defined as T m = [25, 80], α = [1.4,2.5] and η m = [2,14].Fig. 6 shows the region delimited by the cloud of points corresponding to the deviatoric stress invariant σ VM and the volumetric stress invariant I 1 determined for all the elements corresponding to the matrix in the discretized finite element model.Each region in Fig. 6 was obtained according to the respective applied macroscopic failure stress, i.e.X T , X C , Y T , and Y C .It is possible to see in Fig. 6 that the curve obtained for the matrix failure criterion, in terms of the optimal matrix strength parameters, is particularly close to the plotted failure regions.The resulting matrix strength parameters are T m = 43.3MPa,C m = 60.5MPa, and η m = 4.33.

Interface strengths
In order to determine the interfacial strengths, first the interfacial traction along the normal direction t n and the interfacial shear traction t s are determined for all the elements corresponding to the interface in the discretized finite element model.In the present work, all the first elements adjacent to the elements corresponding to the fibers were considered to be members of the interface and their elastic properties were assumed to be equal to the matrix elastic properties.A transformation of stresses is applied over the stresses at microscopic level σ ij in order to obtain the interfacial tractions t i .More details about this transformation can be seen in [40].After computing the interfacial tractions t n and t s for all the elements corresponding to the interface in the discretized finite element model, a cloud of points was plotted on the space t n × t s .Similarly to [44], the interfacial shear strength was assumed as the maximum value of the interfacial shear traction, i.e.
Then, the coordinate of the points located at the boundary of the cloud of points, plotted on the space t n × t s , were determined for all the points whose interfacial traction along the normal direction is positive, i. e. the points with the coordinate t n k > 0. Thus, given Y s and the points , it is possible to obtain Y n according to the Interface Failure criterion defined in Eq. ( 8).According to Eq. ( 8), it is defined as an elliptical relation between t n × t s whose vertex and co-vertex are respectively Y s and Y n .Therefore, depending on the points ( t s k , t n k ) the value of Y ni may vary according to the resulting ellipse.Herein, it is adopted that the interfacial traction strength along the normal direction Y n is obtained using the point which maximizes the co-vertex for the ellipse, i.e.
The assumptions herein adopted for determining the interfacial strengths slightly differ from those presented in [39,44].For instance, Ha et al. [39] mentioned that there is no direct evidence showing the relationship between the interfacial strengths, although they assumed that Y s = 3Y n , with Y n obtained from previous investigations available in the literature.As can be seen in their results, this strategy led to an overestimation of Y s for both of the composite materials investigated in their work.For its part, Macedo et al. [44] considered that Y n = max(t n ) which may eventually lead to an underestimation of Y n .Fig. 7 shows the regions delimited by the cloud of points which correspond to the interfacial tractions along the normal direction t n and to the interfacial shear tractions t s determined for all the elements corresponding to the interface in the discretized finite element model.Analogous to previous plots related to the fiber and matrix strengths, each region in Fig. 7 was obtained according to the respective applied macroscopic failure stress.It is possible to realize from Fig. 7 that the macroscopic longitudinal shear stresses have no effect on the interfacial strengths.Furthermore, it is also possible to verify that the application of macroscopic stresses, which are equivalent to the lamina strengths along

Computation of the micromechanical failure envelopes
Once the strengths of the constituents were determined, the next step involved the application of macroscopic stresses { σ0 } on the RVE.In the present work, it was decided to discretize the stress spaces σ0 ii × σ0 jj and for application of the macroscopic stresses which were defined as where θk = 2πk nS , n S is the total number of directions d Sk
Therefore, given a direction d S k → and an initial point with coordinates , the stresses σ ij at the microscopic level are computed using the Asymptotic Homogenization method for all the elements in the discretized finite element model.Afterwards, whether or not the constituent fails is verified according to the constituent failure criteria described in Section 2. If one or more constituent fails, it is assumed that the material failed and the point with coordinates , is plotted on the material failure envelope for the respective stress space.Otherwise, the macroscopic stresses ) , or ) , is defined.New points for the macroscopic stresses are computed until one or more constituents fail, resulting in ultimate failure of the material.The macroscopic stresses scaling is performed using the Golden Section search technique which is a very reliable line search method (see [67] for more detailed information).For this purpose, an objective function F {σ 0 } is defined and an optimization problem is proposed as with where FFC T is the fiber failure criterion in tension, FFC C is the fiber failure criteria in compression, MFC is the matrix criterion, and IFC is the interface failure criterion defined in Section 2.
Since F {σ 0 } is unimodal within the interval of search, i.e. there is only one minimum point for each direction of search defined on the stress spaces σ0 ii × σ0 jj and σ0 ii × σ0 ij , the convergence is assured and the accuracy of the obtained failure points depends only on the uncertainty defined in the implementation of the method.Once all the directions of search are evaluated, i.e. all the macroscopic stresses on the stress spaces σ0 ii × σ0 jj and σ0 ii × σ0 ij that cause the material failure are determined, the complete failure envelope for the material is obtained.Fig. 8 summarizes the procedure previously described to compute the micromechanical failure envelopes of 3D-printed CFRTPCs.

Special conditions for macroscopic deviatoric stresses in matrix failure
As previously mentioned in Section 4.2, the macroscopic failure stresses equivalent to the longitudinal shear strength S 12 was not taken into account in the determination of the matrix strengths, based on the assumptions for matrix failure described in [39].As a consequence, the failure under pure longitudinal shear load in their work consistently occurred below the longitudinal shear strength obtained experimentally.Furthermore, under combined longitudinal shear and transverse tensile/compressive loads, the failure envelopes were underestimated when compared to the experimental bi-axial testing data of their investigated material.On the other hand, the adoption of the macroscopic failure stresses equivalent to the longitudinal shear strength S 12 for the determination of matrix parameters in [44] resulted in failure envelopes which agreed with the experimental result for the uniaxial ply shear strength.However, it can be verified from their results that the ply failure under transverse compressive loads, with σ 12 = 0, predominantly occurred due to interface failure in view of the overestimated compressive strength.
From the previous discussion, a strategy which accounts for the effects of longitudinal shear and lies between the approaches previously mentioned, should be considered, i.e. the longitudinal shear stresses should be taken into account for the computation of the failure envelopes but the matrix failure criterion should be essentially dominated by transverse tensile and transverse compressive stress components.In this context, the present work aims to propose the adoption of a factor that scales the shear stresses at macroscopic level and defines a relation between transverse tensile stress σ22 and the longitudinal shear stress σ12 for the matrix failure.Assuming a plane stress-state and taking into account the matrix failure criterion curve defined in Section 2.2 (see Eq.
Thus, the ratio F dt1 that results in a ply failure under longitudinal shear stress coincident to the longitudinal shear strength can be written as In principle, the ratio F dt1 applied over the macroscopic longitudinal shear stress solves the problem of the underestimated failure envelope obtained in [39].However, the transition between transverse tension and longitudinal shear can be compromised.In order to circumvent this, an additional interaction factor F dt2 between transverse tensile stresses and longitudinal shear stresses, when σ0 22 > 0, is herein proposed as where β is Taking into account Eqs. ( 25) and ( 26) inside their interval of application, a general interaction expression F dt , in terms of F dt1 and F dt2 can be defined as where sign() is a signum function that returns Therefore, the macroscopic longitudinal shear stress applied to the RVE becomes The adoption of F dt leads to a progressive effect of σ12 on the matrix failure where the maximum effect is reached when σ22 = 0. Consequently, a better transition is expected between tensile and longitudinal shear failure stresses for the matrix in addition to the observance of the uniaxial ply failure strengths obtained from the experimental test data.

Results and discussion
Below the failure envelopes are presented for the 3D-printed CFRTPCs herein investigated.The analytical failure envelopes obtained from early and recent analytical failure theories, i.e.Tsai-Wu [54], Azzi-Tsai [68], Puck and Schürmann [69,70] and Gu-Chen [71], are also plotted.Since in [26] an expansion of the Puck and Schürmann Inter-Fiber Fracture criterion is proposed for 3D-printed materials, this expansion, also referred to as ExPan, is applied for the failure plane σ 22 × τ 12 , instead of the original Puck and Schürmann.Fig. 9 displays the computed failure envelopes on the stress space σ 22 × τ 12 for the 3D-printed CFRTPCs.In a general context, it has been verified in Fig. 9 that the micromechanical failure envelope seems to be in good agreement with the analytical failure envelopes when the material is loaded in tension.Part of this good agreement can be attributed to the adoption of F dt which led to a progressive effect of σ12 on the matrix failure.As expected, when compressive macroscopic stresses are applied, differences between the failure envelopes are verified.However, the failure predicted by the micromechanical approach was mainly   governed by the matrix failure.Fig. 10 displays the computed failure envelopes in the stress space σ 11 × σ 22 for the 3D-printed CFRTPCs.In this stress space, the herein proposed micromechanical approach provides a failure envelope which seems to behave as an "average-intersective approach" since it is slightly placed in the middle of the intersection region between the different analytical approaches.After a close analysis, it can be seen that, under bi-axial loads with longitudinal tensile stresses, the failure is systematically predicted to occur in the matrix.Fig. 11 shows the computed failure envelopes on the stress space σ 11 × τ 12 for the 3D-printed CFRTPCs.It is possible to verify from Fig. 11 that the micromechanical approach provided a failure envelope which is considerably less conservative than those provided by the early failure theories.However, as expected, it is more conservative than the failure envelope predicted by the original Puck and Schürmann theory, since on this plane the latter is highly similar to the maximum stress theory.
In order to provide more elements into the discussion, a geometrical analysis based on the area occupied by the failure envelopes was performed using the micromechanical failure envelopes as the reference.It is worth mentioning that it could not be found in the available literature reliable bi-axial testing data for reinforced and unreinforced 3D-printed materials.Therefore, it is believed that geometrical comparisons between the failure envelopes is an adequate first step to quantify their difference.Table 3 summarizes the area differences.It is possible to verify that the micromechanical failure envelope on the stress space σ 22 × τ 12 is surprisingly less conservative than the other failure envelopes.However, it is considerably close to that provided by the expansion of the Puck and Schürmann proposed in [26].The previously mentioned "average-intersective approach" seen on the stress space σ 11 × σ 22 becomes more evident analyzing the results in Table 3.The quantified differences in Table 3 show that the difference between the failure envelopes are not negligible, mostly when micromechanical failure envelopes are compared to those provided by the early failure theories.These results support the hypothesis that mechanism based techniques should be applied in order to better predict the failure of 3D-printed composite materials.

Conclusion
The present work was aimed at conducting a micromechanical failure analysis of carbon fiber reinforced 3D-printed composite materials.For this purpose, the failure at the constituent level was investigated according to different failure criteria defined for each one of the constituent phases, i.e. fiber, matrix, and interface.It was then assumed that, when one of the previous constituents initially failed, the 3Dprinted composite material was considered to have failed.In regards to the matrix, the present work demonstrated that, considering macroscopic stresses which are equivalent to the lamina longitudinal shear strength in the determination of matrix strengths, leads to an overestimation of the matrix failure strength in compression.As a consequence, when subjected to macroscopic stresses equivalent to the lamina compressive strength, the material failure occurs in the interface.Notwithstanding, it was also demonstrated that, neglecting the lamina longitudinal shear strength in the determination of matrix strengths, the matrix failure strength under longitudinal shear stresses is highly underestimated.In order to circumvent this, the present work proposed a different strategy based on the adoption of a factor that scales the longitudinal shear stresses at macroscopic level.An additional relation between the transverse tensile stress and the longitudinal shear stress was also determined for the matrix failure.The failure envelopes obtained by the micromechanical approach showed its ability to predict the failure of 3D-printed CFRTPCs.Furthermore, it is believed that the characteristics inherent to the microstructure, which are implicitly relayed to the RVE, contributed to obtaining a more realistic failure envelope, differing at specific regions from those obtained by analytical approaches which are not sensitive to fiber arrangement.It could be seen from the available literature that, when the material is subjected to a bi-axial tensile stress-state, the matrix failure may be governed by the dilatational energy.However, preliminary studies have demonstrated that the difference between the results obtained by the stress-based approach and those obtained by the dilatational energy-based approach are practically the same.Additionally, it is worth noting that a failure criterion accounting for a damage initiation on the interface, and its eventual propagation on the matrix, in combination with a model which considers the visco-plasticity effects commonly seen in thermoplastic matrices, could be very suitable in future analyses.Lastly, it can be said that there still exists the need for more experimental work in order to provide reliable bi-axial testing data, which would support improvements on the approaches respectively adopted in the present work.

}
equivalent to the uniaxial strengths of the lamina, i.e.X T , X C , Y T , Y C and S 12 , are selectively applied to the RVE as illustrated in Fig.4(b).This means that each one of the lamina failure stresses in Fig.4(b) corresponds to one load case applied to the

Fig. 2 .
Fig. 2. Comparison of experimental mechanical properties with homogenized elastic properties computed for the respective RVE size.

Fig. 3 .
Fig. 3. Homogenized mechanical properties obtained for different sizes of RVEs: extension elastic moduli (a), shear elastic moduli (b) and Poisson ratios (c).It is worth noting that in (a) the longitudinal modulus, E 1 , is related to the left axis and the transverse moduli, E 2 and E 3 , are related to the right axis.

Fig. 4 .
Fig. 4. Determination of constituents strengths in 3D-printed CFRTPCs: algorithm applied (a) and application of lamina failure stresses, equivalent to uniaxial strengths, on the numerical model (b).

Fig. 5 .
Fig. 5. Tensile (a) and compressive (b) strengths of the continuous fibers reinforcing thermoplastics in the 3D-printed CFRTPCs herein investigated.FFC T is the acronym for Fiber Failure Criterion in tension and FFC C is the acronym for Fiber Failure Criterion in compression.

Fig. 6 .
Fig. 6.Determination of matrix strength parameters C m , T m , and η m according to the adopted matrix failure criterion.MFC is the acronym for Matrix Failure Criterion.

Fig. 7 .
Fig. 7. Interfacial normal and shear traction strengths of the 3D-printed CFRTPCs herein investigated.IFC is the acronym for interface failure criterion.
T.A.Dutra et al.   the fiber direction, do not cause interfacial failure.The resulting interfacial strengths are computed as Y n = 54.6 MPa and Y s = 36.5MPa.

Table 1
Summary of minimum recommended size δ according to the adopted criteria.

Table 3
Comparison between the area of the micromechanical failure envelopes and the area of the analytical failure envelopes on the stress spaces σ 22 × τ 12 , σ 11 × σ 22 and σ 11 × τ 12 .The difference Δ A was computed according to Δ A = 100 × (A MMF − A AF )/A MMF .