Microstructural Analysis of the Transverse and Shear Behavior of Additively Manufactured CFRP Composite RVEs Based on the Phase-Field Fracture Theory

: Due to the versatility of its implementation, additive manufacturing has become the enabling technology in the research and development of innovative engineering components. However, many experimental studies have shown inconsistent results and have highlighted multiple defects in the materials’ structure thus bringing the adoption of the additive manufacturing method in practical engineering applications into question, yet limited work has been carried out in the material modelling of such cases. In order to account for the effects of the accumulated defects, a micromechanical analysis based on the representative volume element has been considered, and phase-ﬁeld modelling has been adopted to model the effects of inter-ﬁber cracking. The 3D models of representative volume elements were developed in the Abaqus environment based on the ﬁber dimensions and content acquired using machine learning algorithms, while fulﬁlling both geometric and material periodicity. Furthermore, the periodic boundary conditions were assumed for each of the representative volume elements in transversal and in-plane shear test cases,. The analysis was conducted by adopting an open-source UMAT subroutine, where the phase-ﬁeld balance equation was related to the readily available heat transfer equation from Abaqus, avoiding the necessity for a dedicated user-deﬁned element thus enabling the adoption of the standard elements and features available in the Abaqus CAE environment. The model was tested on three representative volume element sizes and the interface properties were calibrated according to the experimentally acquired results for continuous carbon-ﬁber-reinforced composites subjected to transverse tensile and shear loads. This investigation conﬁrmed the consistency between the experimental results and the numerical solutions acquired using a phase-ﬁeld fracture approach for the transverse tensile and shear behavior of additively manufactured continuous-ﬁber-reinforced composites, while showing dependence on the representative volume element type for distinctive load cases.


Introduction
Additive manufacturing (AM) technology of polymeric materials has more frequently been adopted in engineering applications [1,2], especially for prototyping in research and development [2,3], with prospects in mold tooling and tool jigs, limited-run production parts, service parts, maintenance and repair, and lightweight solutions in automotive industries [4][5][6], as well as in medical applications [7][8][9].A comprehensive review aiming to synthesize the most relevant studies on carbon-based nano-composites was conducted in [10].The authors discuss the mechanical and electrical properties of materials reinforced with carbonaceous nanofillers, such as graphene nanoplatelets and carbon nanotubes, and their compatibility with the requirements for biomedical applications, emphasizing the potential of the specific electrical, thermal, and mechanical properties of CNTs fillers in multifunctional, hight-strength, lightweight, and self-lubricating wear resistant composites.In addition, a comprehensive review on graphene nanofillers in metal matrix composites for biomedical applications was conducted in [11], highlighting the influences of nano-additive phases on cytotoxicity and the tailored mechanical properties achieved in implants by adopting scaffold structures.All thing considered, the application of AM technology was also further extended by the introduction of particle or fiber fillers to achieve the beneficial properties of these reinforcing materials.The prospects of embedding nanofillers in additively manufactured composite materials were investigated in [12], highlighting the increased strength and stiffness achieved when the components are reinforced with fibers and fillers, which led to a significant amount of research being concentrated on material selection and enhancing the characteristics of cost-effective printed components [12].However, the experimental observations on the behavior of these enhanced materials show inconsistency with the performance of their conventionally manufactured counterparts [13][14][15], which has also been found as one of the major obstacles to the adoption of AM in the current industries [16].A comprehensive review of the application of AM composites has been revied in [17], highlighting its potential in medical applications.A comprehensive review of the application of AM composites has been revied in [17], while the enhancement of the dielectric and thermal properties of polymeric composites by the addition of fillers was investigated in [18,19] respectively.The mechanical properties of carbon-fiber-reinforced specimens containing various LSSs were studied in [20].While these reinforced materials exhibit heterogeneity, the AM components manifest in manufacturing inherent deposition deficiencies which are often attributed to the microscopic voids and inclusion [15,21,22].Among these defects, the AM parts also depend on various controllable manufacturing parameters as well as the mechanical properties of the constituents, which often results in data inconsistencies and issues uncertainty in the prediction of the material behavior.These FDM-induced flaws were further studied in [21][22][23], indicating the necessity for microscopic analysis and NDT component examinations, as well as quantification of the deficiencies which affect the mechanical behavior of AM CFRP composites; however, the modelling approaches which account for the effects of the AM process are limited.
In order to calculate the resulting material properties in novel composite AM extrusion mixtures, microstructural analysis of these heterogenic materials is essential.This can be achieved by analyzing a representative volume element (RVE) designed according to the available microstructural information on the constituents' geometries and behaviors, as shown in [24][25][26][27], as well as being proven valid in multiple studies [28][29][30][31][32][33].The constituents' geometry can be acquired from 2D microscopy [13,34] and then prepared using AI-thresholding [35] for statistical evaluations, based on which, and some assumptions, the RVE is built or it can be analyzed from the 3D µCT scans [36,37] using the voxel approach [38].Such models are prepared in a CAE environment where each of the constituents adopts its constitutive behavior.Since these material models are usually acquired from macroscale analyses, additional assumptions on the intermaterial contact behavior [13,28,39,40] have to be adopted.In commercial FEA software, these interactions are implemented using cohesive elements [28] or alternatively using phase-field modelling [39][40][41].
The phase field approach was initially developed to overcome the problems of moving interface boundaries and topology adjustments when simulating the morphologies of evolving interfaces in cases of merging or division [42,43].Based on the initial formulations presented in [44,45], the model uses a set of field variables Φ which are described using partial differential equations (PDEs) and they distinguish the two phases (0 and 1), between which a smooth transition is assumed when calculated in the proximity of the interface [46].These sets of variables are continuous functions of time and spatial coordinates and can be integrated without explicit treatment of interface conditions [42,43].By adopting the thermodynamic free energy approach, the phase field evolution can be defined in terms of temperature, strain, or concentration to predict various changes in microstructural features, leading to frequent application in simulations of microstructural evolutions [47] and the modelling of fracture mechanics [48].
Furthermore, the development of solution algorithms for the coupled deformationfracture problem was presented in [49,57,78,79].The total potential energy is minimized according to the displacement field u and phase field Φ, while Φ is solved as a damage variable [80] at the finite element node instead of at the integration point as in cases of local damage models [46].This is achieved using a user-defined element (UEL) prior to the implementation of the user defined material model (UMAT), which limits the modelling and visualization capabilities of commercial software [46].This has been accounted for in [46,80] by adopting the similarity between the phase field balance and heat transfer equations which enables the calculation of phase-field damage Φ variable as an additional degree of freedom without the necessity for the UEL mesh definition.To retain all of the the Abaqus built-in modelling and meshing features, this approach has been implemented for the micromechanical damage modelling of AM carbon fiber-reinforced thermoplastic polymer (CFRP) composite RVE and validated on experimentally acquired data.All things considered, this work is focused on the application of the phase field approach in modelling the microstructural failure of AM CFRP composite RVE subjected to tensile and shear loads.Therefore, as an extension to the previously published work where the RVEs were modelled based on continuum mechanics and cohesive finer/matrix interfaces [13,81], the behavior of the AM CFRP composite was modelled in this work using a phase field modelling approach, introducing microscopic damage and showing better scalability of the RVE based on microscopic imaging in contrast to the previously acquired results in [13,81].

Materials and Methods
The consistency of an AM material may vary due to the ratios and microscopic properties of the constituents, the chemical properties of the sizing agents, layer thickness, nozzle diameter, humidity, filling type, and density, as well as the extrusion speed and temperature [82].Since many of these parameters could be controlled, an iterative design process in novel material development is essential.Therefore, a micromechanical modelling approach based on the representative volume element has been proposed in this work.
The CFRP composite specimens for microscopic analysis were designed as 12 layered 25 × 25 mm laminated plates with a cross-ply stacking sequence equal to [0/90 2 /0/90 2 ] s , as shown in Figure 1, and additively manufactured using a Markforged X7 3D printer.This LSS was adopted after a preliminary microstructural examination of similar AM composites.A similar LSS is used in acquiring the multiaxial response through uniaxial tensile tests on rectangular specimens [83], which was also studied in [14].This sequence also enables the measurements of unidirectionally reinforced layer composition and gives a better insight into laminar dimensions, layer-wise defect accumulation, and repeatability issues between single, double, and multiple adjacent layers.However, these specimens did not account for the lengthwise defect accumulation due to uneven material solidification.
composites.A similar LSS is used in acquiring the multiaxial response through uniaxial tensile tests on rectangular specimens [83], which was also studied in [14].This sequence also enables the measurements of unidirectionally reinforced layer composition and gives a better insight into laminar dimensions, layer-wise defect accumulation, and repeatability issues between single, double, and multiple adjacent layers.However, these specimens did not account for the lengthwise defect accumulation due to uneven material solidification.
(a) (b) The specimens were then cut along and perpendicular to the dominant fiber direction in order to acquire x-z and y-z cross-sections, while the x-y cross-section is acquired by grinding the polymeric top layer of the specimens.The cross-sections were polished and prepared for observation using the FEI-QUANTA-250-FEG scanning electron microscope.The acquired images were inspected and, due to the ununiform matrix background, the images were segmented based on pixel classification using Trainable WEKA segmentation machine learning algorithms which are readily available in FIJI [35] opensource software.The acquired results were statistically evaluated and compared with the data available in the relevant literature.Following the comparison, the RVEs were designed according to [13,31,84], while adopting the microscopic properties acquired from the four middle layers.The first case (RVE-1) was designed as a single fiber unit cell; the rectangular and hexagonal fiber arrangements have been adopted for the second (RVE-2) and the third (RVE-3), respectively, while keeping the smallest possible RVE size for the measured fiber volume ratio to increase the computational efficiency.However, to acquire more accurate results from the analysis of small-size RVEs, the fibers were placed at the RVE boundary edges [84], while ensuring geometrical, material, and mesh periodicity [13,31].This was adopted for the cases of rectangular and hexagonal fiber arrays, while for the unit-cell case, both conditions could not be fulfilled.
In order to distinguish between the fiber, the matrix, and the fiber/matrix interface behaviors, each RVE domain was divided into three subdomains.The model was discretized using coupled temperature-displacement hexahedral elements (C3D8T), with the predefined temperature field indicating the initial state of damage within the material.Furthermore, the linear elastic behavior was assumed for the phase field damage framework [33], while the fiber and matrix's mechanical properties were adopted from [85] and [86], respectively; as shown in Table 1.The constituents' properties were adopted from the literature since the commercially available data are available only for the UD composite reinforced in the direction of the layers with the elastic modulus of 60 GPa and tensile strength of 800 MPa [87], while the matrix data are available for the AM polymer used in the secondary nozzle, and not in the analyzed CFRP material.The toughness K1C was adopted for polyamide according to [88] and recalculated to acquire the GC values.In The specimens were then cut along and perpendicular to the dominant fiber direction in order to acquire x-z and y-z cross-sections, while the x-y cross-section is acquired by grinding the polymeric top layer of the specimens.The cross-sections were polished and prepared for observation using the FEI-QUANTA-250-FEG scanning electron microscope.The acquired images were inspected and, due to the ununiform matrix background, the images were segmented based on pixel classification using Trainable WEKA segmentation machine learning algorithms which are readily available in FIJI [35] open-source software.The acquired results were statistically evaluated and compared with the data available in the relevant literature.Following the comparison, the RVEs were designed according to [13,31,84], while adopting the microscopic properties acquired from the four middle layers.The first case (RVE-1) was designed as a single fiber unit cell; the rectangular and hexagonal fiber arrangements have been adopted for the second (RVE-2) and the third (RVE-3), respectively, while keeping the smallest possible RVE size for the measured fiber volume ratio to increase the computational efficiency.However, to acquire more accurate results from the analysis of small-size RVEs, the fibers were placed at the RVE boundary edges [84], while ensuring geometrical, material, and mesh periodicity [13,31].This was adopted for the cases of rectangular and hexagonal fiber arrays, while for the unit-cell case, both conditions could not be fulfilled.
In order to distinguish between the fiber, the matrix, and the fiber/matrix interface behaviors, each RVE domain was divided into three subdomains.The model was discretized using coupled temperature-displacement hexahedral elements (C3D8T), with the predefined temperature field indicating the initial state of damage within the material.Furthermore, the linear elastic behavior was assumed for the phase field damage framework [33], while the fiber and matrix's mechanical properties were adopted from [85] and [86], respectively; as shown in Table 1.The constituents' properties were adopted from the literature since the commercially available data are available only for the UD composite reinforced in the direction of the layers with the elastic modulus of 60 GPa and tensile strength of 800 MPa [87], while the matrix data are available for the AM polymer used in the secondary nozzle, and not in the analyzed CFRP material.The toughness K 1C was adopted for polyamide according to [88] and recalculated to acquire the G C values.In contrast to the fiber and the matrix, the values for the interface properties were assumed and calibrated according to the experimentally observed material behavior.
In order to appropriately simulate the behavior of the surrounding material, periodic boundary conditions (PBC) have been enforced on the RVEs boundary corners, edges, and surfaces [28,[30][31][32][33], where the PBCs and the node constrain equations were defined using the automated procedure presented in [31,32].However, instead of focusing on extracting the elastic properties, the procedure was modified to calculate the material degradation and damage using user-defined subroutines and simulate the complete material response until failure.
The mesh sensitivity was studied on three distinctive cases of RVEs with assumed hexahedral fiber placement.The curvature maximal deviation factor h/L was assumed as 0.01 for each of the tested cases, while the approximate element sizes were adopted as 2.5 × 10 −4 , 3 × 10 −4 , 4 × 10 −4 , and 5 × 10 −4 which resulted in 1,115,403, 60,860, 53,876, and 33,781 elements, respectively.Returning consistent results within the elastic region, the models were discretized using coupled temperature-displacement hexahedral elements (C3D8T) assuming the approximate element size equal to 4 × 10 −4 , resulting in mesh sizes of 23,163, 57,300, and 53,876 elements for RVE-1, RVE-2, and RVE-3, respectively.Finally, the models were subjected to transversal and shear loads and analyzed using a phase field model.
In a standard phase field model formulation (AT2) presented in Equation ( 1), the damage is described using the phase field variable Φ which takes the values between zero for the undamaged and one for the fully damaged material, while the damage evolution depends on the balance between the energy stored within the material and the energy released during fracture [46,80].
In this expression, the fracture toughness of the material and the strain energy density are represented by the variables G c , and ψ, respectively, while the variable l stands for the phase field length which is used to define the size of the damaged region.In Abaqus, this differential equation is solved by introducing the user-defined element (UEL) before the material model subroutine is called, thus losing the Abaqus built-in modelling, meshing, and visualization features in the process.To account for this, a solution based on the similarities between the phase field damage evolution and the heat transfer problem under steady-state conditions has been presented in [46,80].In this framework, the change in the temperature T for a material with thermal conductivity k which is subjected to a heat source r is calculated according to Equation ( 2), and since this scheme is readily available as an Abaqus built-in feature, it was utilized for solving the differential equation for the phase field damage evolution [46,80].
The complete mathematical formulation of the phase field theory is presented in [46,80,89] as well as in Appendix A, and within this framework, the initial temperature is given as a predefined field variable, the material thermal conductivity is equal to one, while the variable r is redefined according to the phase field fracture theory and introduced using the UMAT subroutine.Following these assumptions, the initial temperature is usually given equal to zero and it describes the initial state of damage, which evolves according to the phase field damage formulation presented as r until the fully damaged state is reached.Within the presented study, the constituents were adopted as isotropic, while presuming linear elastic behavior based on four mechanical properties for each of the constituents including the modulus of elasticity E, the Poisson ratio ν, the characteristic phase field length scale l, and the fracture toughness G c .Besides these variables, the model variations are also presented in [46,80] including AT1 [73], AT2 [45], and the phase field cohesive approach PF-CZM [74], which also takes the material tensile strength into account.Furthermore, the monolithic and staggered solution schemes are also presented in [46,80], as well as a hybrid [90] or anisotropic [78] approach for the volumetric-deviatoric [91], and the spectral [78] strain energy split schemes.Within the framework of the presented research, the standard AT2 model has been adopted for each of the constituents in all of the RVE cases.While using the staggered solution scheme without strain decomposition in order to achieve convergence with a larger number of iterations, the Abaqus general solution controls I 0 , I R , I P , I C , I L , and I G have been set to 5000 [46,80].The RVEs were analyzed using the adopted UMAT subroutine to be compared with the experimental results.
Therefore, two sets of continuous fiber-reinforced specimens with a polyamide matrix have been designed as flat rectangular laminates, as shown in Figure 2.
J. Compos.Sci.2023, 7, x FOR PEER REVIEW 7 of 24 damage formulation presented as r until the fully damaged state is reached.Within the presented study, the constituents were adopted as isotropic, while presuming linear elastic behavior based on four mechanical properties for each of the constituents including the modulus of elasticity , the Poisson ratio , the characteristic phase field length scale , and the fracture toughness  .Besides these variables, the model variations are also presented in [46,80] including AT1 [73], AT2 [45], and the phase field cohesive approach PF-CZM [74], which also takes the material tensile strength into account.Furthermore, the monolithic and staggered solution schemes are also presented in [46,80], as well as a hybrid [90] or anisotropic [78] approach for the volumetric-deviatoric [91], and the spectral [78] strain energy split schemes.Within the framework of the presented research, the standard AT2 model has been adopted for each of the constituents in all of the RVE cases.While using the staggered solution scheme without strain decomposition in order to achieve convergence with a larger number of iterations, the Abaqus general solution controls  ,  ,  ,  ,  , and  have been set to 5000 [46,80].The RVEs were analyzed using the adopted UMAT subroutine to be compared with the experimental results.Therefore, two sets of continuous fiber-reinforced specimens with a polyamide matrix have been designed as flat rectangular laminates, as shown in Figure 2.  The specimens in the first set UD-90 were designed according to the ASTM-D3039 [92] standard as transversely reinforced, unidirectionally reinforced laminates with an LSS equal to [90] , while for the second set SH-45, according to the ASTM-D3518 [93] standard, they were designed as multidirectionally reinforced laminates with an LSS equal to [45/−45] which resolves the uniaxially applied tensile loads into in-plane shear loads; the visual representation of the adopted LSS is for each specimen set presented in Figure 3.The specimens in the first set UD-90 were designed according to the ASTM-D3039 [92] standard as transversely reinforced, unidirectionally reinforced laminates with an LSS equal to [90] 8 , while for the second set SH-45, according to the ASTM-D3518 [93] standard, they were designed as multidirectionally reinforced laminates with an LSS equal to [45/ − 45] 4S which resolves the uniaxially applied tensile loads into in-plane shear loads; the visual representation of the adopted LSS is for each specimen set presented in Figure 3.Both sets of specimens were tested uniaxially in quasistatic conditions with v = 0.01 mm/s using a hydraulic testing system and monitored with both contact and optical extensometers.The measurements have been taken on the opposite sides of the specimens, monitoring the strains of one surface using the epsilon-tech axial extensometer Model-3542, and the opposite one with a GOM Aramis 5M (GigE) adjustable base system with 35 mm lenses.The image resolution was 2448 × 2050 pixels, while the cameras were positioned at 560 mm from the specimen, with the distance between the cameras being equal to 265 mm, closing the angle of 26°.The sensor calibration was conducted using the software-guided protocol at 22.5 °C resulting in an average deviation of 0.048 pixels in a measuring volume of 130 × 110 mm × 90 mm.Such a measurement protocol was selected to check the consistency of the measurements through the specimen thickness and to detect excessive delamination within the gauge length.The measurements have been systematized for comparison with the RVE results and they are presented in the following chapter.

Results
The results are given as the microstructural observations and macroscopic experimental results, following by design of the RVEs and numerical analysis.

Microstructural Inspection
In order to obtain the necessary input for the RVE design, the material microstructure was examined on multiple cross-sections in the y-z, x-z, and x-y planes.The cross-sections were polished and scanned using an FEI-QUANTA-250-FEG scanning electron microscope; the results are presented in Figure 4.Both sets of specimens were tested uniaxially in quasistatic conditions with v = 0.01 mm/s using a hydraulic testing system and monitored with both contact and optical extensometers.The measurements have been taken on the opposite sides of the specimens, monitoring the strains of one surface using the epsilon-tech axial extensometer Model-3542, and the opposite one with a GOM Aramis 5M (GigE) adjustable base system with 35 mm lenses.The image resolution was 2448 × 2050 pixels, while the cameras were positioned at 560 mm from the specimen, with the distance between the cameras being equal to 265 mm, closing the angle of 26 • .The sensor calibration was conducted using the software-guided protocol at 22.5 • C resulting in an average deviation of 0.048 pixels in a measuring volume of 130 × 110 mm × 90 mm.Such a measurement protocol was selected to check the consistency of the measurements through the specimen thickness and to detect excessive delamination within the gauge length.The measurements have been systematized for comparison with the RVE results and they are presented in the following chapter.

Results
The results are given as the microstructural observations and macroscopic experimental results, following by design of the RVEs and numerical analysis.

Microstructural Inspection
In order to obtain the necessary input for the RVE design, the material microstructure was examined on multiple cross-sections in the y-z, x-z, and x-y planes.The cross-sections were polished and scanned using an FEI-QUANTA-250-FEG scanning electron microscope; the results are presented in Figure 4.Both sets of specimens were tested uniaxially in quasistatic conditions with v = 0.01 mm/s using a hydraulic testing system and monitored with both contact and optical extensometers.The measurements have been taken on the opposite sides of the specimens, monitoring the strains of one surface using the epsilon-tech axial extensometer Model-3542, and the opposite one with a GOM Aramis 5M (GigE) adjustable base system with 35 mm lenses.The image resolution was 2448 × 2050 pixels, while the cameras were positioned at 560 mm from the specimen, with the distance between the cameras being equal to 265 mm, closing the angle of 26°.The sensor calibration was conducted using the software-guided protocol at 22.5 °C resulting in an average deviation of 0.048 pixels in a measuring volume of 130 × 110 mm × 90 mm.Such a measurement protocol was selected to check the consistency of the measurements through the specimen thickness and to detect excessive delamination within the gauge length.The measurements have been systematized for comparison with the RVE results and they are presented in the following chapter.

Results
The results are given as the microstructural observations and macroscopic experimental results, following by design of the RVEs and numerical analysis.

Microstructural Inspection
In order to obtain the necessary input for the RVE design, the material microstructure was examined on multiple cross-sections in the y-z, x-z, and x-y planes.The cross-sections were polished and scanned using an FEI-QUANTA-250-FEG scanning electron microscope; the results are presented in Figure 4.The acquired images were inspected and, due to the ununiform matrix background, the images were segmented based on pixel classification using Trainable WEKA segmentation machine learning algorithms in FIJI as shown in Figure 5.In the classification process, three sets of categories have been defined: the fibers, matrix, and voids with debris, each associated with the appropriate pixels.This procedure was conducted on ten randomly selected areas, while the number of included layers depended on the particular crosssection.
The acquired images were inspected and, due to the ununiform matrix backgro the images were segmented based on pixel classification using Trainable W segmentation machine learning algorithms in FIJI as shown in Figure 5.In classification process, three sets of categories have been defined: the fibers, matrix, voids with debris, each associated with the appropriate pixels.This procedure conducted on ten randomly selected areas, while the number of included layers depen on the particular cross-section.The Figure 5b shows the classified microstructure after being analyzed based on p recognition using the machine learning algorithms.The coloured regions define detected phases after the training sessions and correspond to the constituents: the fi matrix, and voids and debris present within the analyzed material.Furthermore constituent's volume ratios were identified using the colour threshold adjustmen multiple randomly selected images from each of the analyzed cross-sections, where major contribution to the result was the images containing four middle layers in the cross-section, as shown in Figure 5.Following the classification algorithm, the fibers the ratios of the fibers, matrix, and debris and voids were identified for each of the sele images from which the average value was acquired and is presented in Table 2. Howe due to the irregularity of the fiber arrangement, the local fiber volume ratio was meas on ten randomly selected 100 µm × 100 µm zones within the middle four-layer showing consistency with the measured average, while being inconsistent with constituents' volume ratios found in material agglomerations.In contrast to the vol ratios, the geometrical variable such as the fiber diameter and layer height could been measured directly from the acquired images, and additional image refinements w not necessary for those datasets.As presented in Figure 6, the statistical analysis o fiber diameters showed consistency with the normal distribution, while larger discrepancies from the normal distribution could be observed in the layer height distr tion which is also indicated by the Anderson-Darling value lower than 0.005.The Figure 5b shows the classified microstructure after being analyzed based on pixel recognition using the machine learning algorithms.The coloured regions define the detected phases after the training sessions and correspond to the constituents: the fibers, matrix, and voids and debris present within the analyzed material.Furthermore, the constituent's volume ratios were identified using the colour threshold adjustment on multiple randomly selected images from each of the analyzed cross-sections, where the major contribution to the result was the images containing four middle layers in the y-z cross-section, as shown in Figure 5.Following the classification algorithm, the fibers and the ratios of the fibers, matrix, and debris and voids were identified for each of the selected images from which the average value was acquired and is presented in Table 2.However, due to the irregularity of the fiber arrangement, the local fiber volume ratio was measured on ten randomly selected 100 µm × 100 µm zones within the middle four-layer area, showing consistency with the measured average, while being inconsistent with the constituents' volume ratios found in material agglomerations.In contrast to the volume ratios, the geometrical variable such as the fiber diameter and layer height could have been measured directly from the acquired images, and additional image refinements were not necessary for those datasets.As presented in Figure 6, the statistical analysis of the fiber diameters showed consistency with the normal distribution, while larger data discrepancies from the normal distribution could be observed in the layer height distribution which is also indicated by the Anderson-Darling value lower than 0.005.The results of the statistical evaluations have been compared with the v published in the relevant literature, and the differences between the volume ratio acq from the sample using microscopy and from the filament using pyrolysis are highlig in Table 2.The comparison yielded consistency for the values of fiber diameters, w significant discrepancies were documented for the total and local fiber volume ratios different machines and parameters were used in manufacturing.Conclusively, despite the microstructural inspection having been proven vali the identification of the constituents' properties, it also disclosed the inconsi repeatability of layer thickness in the AM process, which is seldom discussed essential for the performance of fiber-reinforced composite laminates.Despite comparison with the relevant data in the literature yielding consistency for the valu the fiber diameters, significant discrepancies were observed for the total and local volume ratios.These discrepancies could be influenced by the distinctive machines parameters used in the manufacturing process but they were also expected sinc microstructural images were extracted from the laminated specimen [13], instead o filament presented in [85].

Experimental Acquisition of Lamina Properties
In order to identify the mechanical properties of AM continuous-fiber-reinfo composites, experimental studies on transversal and shear behavior have been condu following the guidelines presented in ASTM-D3039 [92] and ASTM-D3518 [93] stand The dimensions of each specimen were assessed in multiple cross-sections withi The results of the statistical evaluations have been compared with the values published in the relevant literature, and the differences between the volume ratio acquired from the sample using microscopy and from the filament using pyrolysis are highlighted in Table 2.The comparison yielded consistency for the values of fiber diameters, while significant discrepancies were documented for the total and local fiber volume ratios since different machines and parameters were used in manufacturing.
Conclusively, despite the microstructural inspection having been proven valid for the identification of the constituents' properties, it also disclosed the inconsistent repeatability of layer thickness in the AM process, which is seldom discussed, but essential for the performance of fiber-reinforced composite laminates.Despite the comparison with the relevant data in the literature yielding consistency for the values of the fiber diameters, significant discrepancies were observed for the total and local fiber volume ratios.These discrepancies could be influenced by the distinctive machines and parameters used in the manufacturing process but they were also expected since the microstructural images were extracted from the laminated specimen [13], instead of the filament presented in [85].

Experimental Acquisition of Lamina Properties
In order to identify the mechanical properties of AM continuous-fiber-reinforced composites, experimental studies on transversal and shear behavior have been conducted, following the guidelines presented in ASTM-D3039 [92] and ASTM-D3518 [93] standards.The dimensions of each specimen were assessed in multiple cross-sections within the gauge length using a digital micrometer, and the acquired values for the widths and thicknesses were averaged and are reported in Table 3.The acquired values were analyzed, showing inconsistencies in comparison with the CAE data, as presented in Figure 2, caused by the removal of the support material.Therefore, to represent the manufactured material more accurately, the measured values were adopted in the numerical analysis.The specimens were subjected to uniaxial tensile loads and monitored with both contact and non-contact strain measurement techniques.During the experimental protocol, it was confirmed that the behavior of AM UD-90 specimens is highly influenced by the matrix response and also the fact that these specimens exhibit inferior mechanical properties in comparison to their injection-moulded counterparts [86].The material bonding defects influenced by the fiber/matrix interface and the raster-induced defects of the additive manufacturing process act as stress concentrators which can be observed in Figure 7 as localized strains between the deposition paths.Consequently, these effects diminish the strength and stiffness of the AM materials, and due to their stochastic nature, they cause material inconsistencies, which is also confirmed by the scatter in the experimental data, and thus had to be considered in material modelling.The specimens were subjected to uniaxial tensile loads and monitored with both contact and non-contact strain measurement techniques.During the experimental protocol, it was confirmed that the behavior of AM UD-90 specimens is highly influenced by the matrix response and also the fact that these specimens exhibit inferior mechanical properties in comparison to their injection-moulded counterparts [86].The material bonding defects influenced by the fiber/matrix interface and the raster-induced defects of the additive manufacturing process act as stress concentrators which can be observed in Figure 7 as localized strains between the deposition paths.Consequently, these effects diminish the strength and stiffness of the AM materials, and due to their stochastic nature, they cause material inconsistencies, which is also confirmed by the scatter in the experimental data, and thus had to be considered in material modelling.Additionally, the SH-45 specimens represent a specific case of multidirectionally reinforced laminates where the adopted LSS induces multiaxial in-plane shear stress during uniaxial tension, while, similar to the case of UD-90, multiple strain localizations can be observed, as shown in Figure 7.In contrast to the UD-90 case, the SH-45 specimens demonstrated more consistent results up to 2% of shear strain, which is followed by an increase in data scatter up to 5% of shear strain.The anticipated ductile behavior for the SH-45 specimens resolves in the fiber rearrangement phenomenon which can lead to 25% of strain before failure.Since these large strains are not applicable in composite structures, the ASTM-D3518 standard [93] recommends limiting the observations of the experimental results at 5% of shear strain, in case the failure does not occur before.The shear strength in SH-45 specimens was also acquired according to the ASTM standard as the intersection between the shear stress-strain curve and the 0.2% offset of its shear modulus.The experimentally acquired data were systematized and prepared for comparison with the numerical results for each of the RVEs.

RVE Design
According to the acquired microscopic constituents' properties, three sets of RVE models have been developed as presented in Figure 8, where RVE-1 represents a unit cell Additionally, the SH-45 specimens represent a specific case of multidirectionally reinforced laminates where the adopted LSS induces multiaxial in-plane shear stress during uniaxial tension, while, similar to the case of UD-90, multiple strain localizations can be observed, as shown in Figure 7.In contrast to the UD-90 case, the SH-45 specimens demonstrated more consistent results up to 2% of shear strain, which is followed by an increase in data scatter up to 5% of shear strain.The anticipated ductile behavior for the SH-45 specimens resolves in the fiber rearrangement phenomenon which can lead to 25% of strain before failure.Since these large strains are not applicable in composite structures, the ASTM-D3518 standard [93] recommends limiting the observations of the experimental results at 5% of shear strain, in case the failure does not occur before.The shear strength in SH-45 specimens was also acquired according to the ASTM standard as the intersection between the shear stress-strain curve and the 0.2% offset of its shear modulus.The experimentally acquired data were systematized and prepared for comparison with the numerical results for each of the RVEs.

RVE Design
According to the acquired microscopic constituents' properties, three sets of RVE models have been developed as presented in Figure 8, where RVE-1 represents a unit cell containing only one central fiber, while both RVE-2 and RVE-3 contain another fiber which is divided into four quarters placed in the RVE corners to ensure better model predictions using a smaller RVE.Furthermore, to account for the material bonding deficiencies, both the fiber/matrix and the raster-induced bonding inconsistencies have to be accounted for.However, since the scale of raster-induced defects exceeded the minimal RVE dimensions, both the matrix/fiber and the deposition contact deficiencies have been introduced together at the fiber/matrix interface as a distinctive subdomain.
J. Compos.Sci.2023, 7, x FOR PEER REVIEW 13 of containing only one central fiber, while both RVE-2 and RVE-3 contain another fib which is divided into four quarters placed in the RVE corners to ensure better mod predictions using a smaller RVE.Furthermore, to account for the material bondi deficiencies, both the fiber/matrix and the raster-induced bonding inconsistencies have be accounted for.However, since the scale of raster-induced defects exceeded the minim RVE dimensions, both the matrix/fiber and the deposition contact deficiencies have be introduced together at the fiber/matrix interface as a distinctive subdomain.In order to check the significance of the discretization on the numerical results mesh sensitivity analysis was conducted on the RVE-3 case.The models were mesh using hexahedral elements (C3D8T) with characteristic lengths of 2.5 × 10 , 3 × 10 4 × 10 , and 5 × 10 , resulting in 1,115,403, 60,860, 53,876, and 33,781 elemen respectively, returning consistent results within the elastic region, while deviating fro the average values for the post-yielding region in the case of 2.5 × 10 , as presented Figure 9. Since the presented discrepancies were less significant than the scatter within t experimental dataset, the characteristic element length of 4 × 10 was adopted for all the RVE cases.With the mesh size adopted, and the domains discretized, an automatiz procedure [31] was adopted for enforcing the periodic boundary conditions for transve and in-plane shear cases.In order to include the failure analysis of heterogenic materia the procedure was modified to use the phase-field model subroutines presented in [46,8 The necessary mechanical properties for the fiber and matrix materials were adopted fro In order to check the significance of the discretization on the numerical results, a mesh sensitivity analysis was conducted on the RVE-3 case.The models were meshed using hexahedral elements (C3D8T) with characteristic lengths of 2.5 × 10 −4 , 3 × 10 −4 , 4 × 10 −4 , and 5 × 10 −4 , resulting in 1,115,403, 60,860, 53,876, and 33,781 elements, respectively, returning consistent results within the elastic region, while deviating from the average values for the post-yielding region in the case of 2.5 × 10 −4 , as presented in Figure 9.
containing only one central fiber, while both RVE-2 and RVE-3 contain another fiber which is divided into four quarters placed in the RVE corners to ensure better model predictions using a smaller RVE.Furthermore, to account for the material bonding deficiencies, both the fiber/matrix and the raster-induced bonding inconsistencies have to be accounted for.However, since the scale of raster-induced defects exceeded the minimal RVE dimensions, both the matrix/fiber and the deposition contact deficiencies have been introduced together at the fiber/matrix interface as a distinctive subdomain.In order to check the significance of the discretization on the numerical results, a mesh sensitivity analysis was conducted on the RVE-3 case.The models were meshed using hexahedral elements (C3D8T) with characteristic lengths of 2.5 × 10 , 3 × 10 , 4 × 10 , and 5 × 10 , resulting in 1,115,403, 60,860, 53,876, and 33,781 elements, respectively, returning consistent results within the elastic region, while deviating from the average values for the post-yielding region in the case of 2.5 × 10 , as presented in Figure 9. Since the presented discrepancies were less significant than the scatter within the experimental dataset, the characteristic element length of 4 × 10 was adopted for all of the RVE cases.With the mesh size adopted, and the domains discretized, an automatized procedure [31] was adopted for enforcing the periodic boundary conditions for transverse and in-plane shear cases.In order to include the failure analysis of heterogenic materials, the procedure was modified to use the phase-field model subroutines presented in [46,80].The necessary mechanical properties for the fiber and matrix materials were adopted from Since the presented discrepancies were less significant than the scatter within the experimental dataset, the characteristic element length of 4 × 10 −4 was adopted for all of the RVE cases.With the mesh size adopted, and the domains discretized, an automatized procedure [31] was adopted for enforcing the periodic boundary conditions for transverse and in-plane shear cases.In order to include the failure analysis of heterogenic materials, the procedure was modified to use the phase-field model subroutines presented in [46,80].The necessary mechanical properties for the fiber and matrix materials were adopted from the literature, while the properties for the interface were fitted according to the experimental results.The RVEs were analyzed in the ABAQUS CAE environment, with them reaching convergence in each tested case, and the results are compared with the experimentally acquired data as presented in Figure 10, thus showing consistency for the transverse case and deviations from the experiments in the in-plane shear.
the literature, while the properties for the interface were fitted according to the experimental results.The RVEs were analyzed in the ABAQUS CAE environment, with them reaching convergence in each tested case, and the results are compared with the experimentally acquired data as presented in Figure 10, thus showing consistency for the transverse case and deviations from the experiments in the in-plane shear.According to the data presented in Figure 10, there is a consistency between the modelled behaviors for the transverse tensile cases regardless of the RVE type.Furthermore, the model predictions correspond to the experimentally acquired data returning values within the scatter range.However, when subjected to in-plane shear loads, the RVE-1 and RVE-2 diverge from both the experimental results and the results acquired for the RVE-3, underestimating the shear strength and modulus, which in contrast are overestimated by the RVE-3 response.Furthermore, the resulting contour plots of the damage index are presented in Figure 11, thus supporting the interfiber damage initiation.According to the data presented in Figure 10, there is a consistency between the modelled behaviors for the transverse tensile cases regardless of the RVE type.Furthermore, the model predictions correspond to the experimentally acquired data returning values within the scatter range.However, when subjected to in-plane shear loads, the RVE-1 and RVE-2 diverge from both the experimental results and the results acquired for the RVE-3, underestimating the shear strength and modulus, which in contrast are overestimated by the RVE-3 response.Furthermore, the resulting contour plots of the damage index are presented in Figure 11, thus supporting the interfiber damage initiation.
the literature, while the properties for the interface were fitted according to the experimental results.The RVEs were analyzed in the ABAQUS CAE environment, with them reaching convergence in each tested case, and the results are compared with the experimentally acquired data as presented in Figure 10 According to the data presented in Figure 10, there is a consistency between the modelled behaviors for the transverse tensile cases regardless of the RVE type.Furthermore, the model predictions correspond to the experimentally acquired data returning values within the scatter range.However, when subjected to in-plane shear loads, the RVE-1 and RVE-2 diverge from both the experimental results and the results acquired for the RVE-3, underestimating the shear strength and modulus, which in contrast are overestimated by the RVE-3 response.Furthermore, the resulting contour plots of the damage index are presented in Figure 11, thus supporting the interfiber damage initiation.As shown in Figure 11, both RVEs subjected to transverse tensile and shear loads manifest a damage initiation at the fiber/matrix interface, which is then propagated in the matrix region, resulting in fracture.Since the fibers within the RVEs are evenly distributed by the assumption of the hexagonal array, the crack path is significantly influenced by their placements.Therefore, to investigate the crack path evolution, a randomly distributed fiber array should be adopted in further studies.All things considered, it has to be emphasized that despite the transverse micromechanical behavior being well represented, the in-plane shear cases do not follow the experimental results after 1% of shear strain and predict a complete loss of bearing capabilities after exceeding the shear strength calculated according to the ASTM guidelines, which is not the case in these laminates.Therefore, different techniques for the experimental identification and modelling of shear behavior in AM-fiberreinforced composite materials should be considered.

Discussion
In order to predict the behavior of additively manufactured composite materials based on the material microstructural parameters and the constituents' arrangements, a micromechanical modelling approach based on the representative volume element has been proposed in this work.The microscopic analysis was conducted on the representative (25 × 25 mm) cross-ply specimens through multiple cross-sections.Despite acquiring a lot of information on the material microstructure which was validated with the relevant data from the literature and then implemented in the RVE design, it was not possible to recognize any of the deposition paths in either of the cross-sections.This could be caused by the fact that the small specimen size, which is manufactured faster than the previously deposited layer, could have been completely cooled, resulting in a more homogenous microstructure without the characteristic triangular void patterns found in large-scale prints.Therefore, a scale-oriented microstructural inspection should be considered in further studies to address this issue.Furthermore, the 2D microscopic inspection in multiple cross-sections did not offer a complete insight into the materials' microstructure.The reference coordinate system, the printing direction, and misalignments were difficult to establish without assumptions.Therefore, a µCT approach should be considered in the future both for the microstructure and the inspection of failure mechanisms.
In this study, the RVEs were designed assuming various geometrical fiber arrangements with consideration given to the measured fiber volume ratio, while keeping the RVE surfaces symmetrical to enable periodic mesh generation and the implementation of the periodic boundary conditions.Such RVEs fulfill both the geometrical and material periodicity conditions; hence, the fibers complete each other if RVEs' segments are arranged together.However, the geometrical arrangement in these RVEs does not reflect the proper nature of the fiber arrangement found in these materials.Therefore, a statistically random fiber distribution or an input from the µCT scan should be considered in future studies.Three distinctive cases of RVEs were developed in the Abaqus CAE environment based on the results from microscopic inspections.The analysis was conducted using the phase-field fracture model available as an open-source UMAT subroutine, where the phase-field balance equation was related to the heat transfer equation readily available in the Abaqus CAE environment, hence avoiding the necessity for a dedicated UEL.This ensured the adoption of standard elements and features available in Abaqus and therefore the application of automatized procedures for enforcing the periodic boundary conditions for transversal and in-plane shear cases.In this analysis, the linear elastic constituents' behavior was assumed, while the fracture properties were calculated according to the available properties for similar materials.Thus, for better insight into the microstructural behavior, these properties should be measured.
In order to acquire the transversal tensile and in-plane shear mechanical properties of AM CFRP composites, macroscopic evaluations were conducted on unidirectional, transversely reinforced UD-90 and multidirectionally reinforced SH-45 laminates, respectively.Both cases were tested uniaxially, where for SH-45, the adopted LSS resolved the uniaxial tensile stresses into in-plane shear.As anticipated, the transverse behavior was confirmed to be dominated by the matrix's mechanical properties; however, the overall material response was inferior to that of the injection-moulded counterparts found in the literature, thus it was attributed to the fiber/matrix interface and the raster-induced defects of the additive manufacturing process.These effects act as stress concentrators and reduce the strength and stiffness of the AM materials, which was also confirmed in the scattering of the experimentally acquired data.These effects were also documented in the SH-45 specimens, which showed 2% divergence from the average before reaching 5% of shear strain.Therefore, both the fiber/matrix and the raster-induced bonding inconsistencies had to be accounted for in the material modelling.However, as the scale of raster-induced defects exceeded the size of the adopted RVEs, both the matrix/fiber and deposition contact deficiencies were introduced uniformly at the fiber/matrix interface.Since this approach is not physically consistent, yet based on an assumption, a multiscale representation to distinguish these effects should be considered in further studies.However, the properties of the fiber/matrix interface or the material deposition contact need to be acquired experimentally.

Conclusions
Numerical calculations of the transverse tensile and in-plane shear micromechanical behavior of AM CFRP materials based on the material microscopic analysis and the macroscopic response were presented in this study.The procedures used in material manufacturing and preparation for microscopic inspections were discussed, the positive and the negative aspects of the adopted approach were presented, and the key features for acquiring better results in further studies have been highlighted.The RVE design based on geometrical assumptions and supported by the microstructural measurements was also discussed, where many benefits of using µCT over microscopy have been presented.Furthermore, the presented microstructural analysis is based on the proposed analogy between phase-field fracture and the thermal conductivity differential equations since it supports the ABAQUS built-in features, thus enabling the development and analysis of the periodic RVEs using the phase-field fracture theory.This approach simplified the domain divisions based on the various materials and their interfaces, as well as the load application and finally the visualization of the acquired result.Furthermore, both the data scatter and the DIC images confirmed that the inconsistencies in AM CFRP composites are caused in the contact zones of the deposited material.However, since the scale of these deficiencies exceeded the size of the proposed RVEs, the effects were introduced together with the fiber/matrix influence at the constituents' interface.The negative aspects of this approach were discussed, and a multiscale approach with experimentally acquired interface properties was proposed for further study.
All things considered, this study confirmed consistency between the experimental results and the numerical prediction using the phase-field fracture approach for the transverse tensile behavior of AM CFRP composites.The numerical results were independent of the RVE type, returning values within the range of the experimental scatter.However, when analyzing the in-plane shear behavior, RVE-1 and RVE-2 revealed divergence from both the experimental results and the results acquired for the RVE-3, thus underestimating the shear strength and modulus of the tested materials.In contrast, RVE-3 overestimated the mechanical properties in the case of in-plane shear.It was shown that despite the transverse micromechanical behavior being well characterized using this phase-field modelling approach, the in-plane shear cases do not follow the experimental results after 1% of shear strain.Furthermore, after exceeding the shear strength calculated according to the ASTM guidelines, the model predicts a complete specimen failure which does not occur in the experiments.Therefore, different techniques for experimental identification and modelling of shear behavior in AM-fiber-reinforced composite materials should be considered in future work.
Author Contributions: M.G.: writing-original draft, visualization, software, methodology, investigation, and formal analysis.D.L.: writing-review and editing, supervision, project administration, and funding acquisition.M.F.: writing-review and editing, supervision, project administration, and

Figure 3 .
Figure 3. (a) Layer stacking sequence for the UD-90 set, and (b) layer stacking sequence for the SH-45 set; recreated according to [81].

Figure 3 .
Figure 3. (a) Layer stacking sequence for the UD-90 set, and (b) layer stacking sequence for the SH-45 set; recreated according to [81].

Figure 5 .
Figure 5. CFRP laminate's four middle unidirectional layers: (a) SEM image, (b) W classification, (c) the fiber vs. matrix probability map, and (d) voids and debris probability ma

Figure 9 .
Figure 9.Comparison between the experimental results and the RVE-3 response for characteri element length in the range from 2.5 × 10 to 5 × 10 .

Figure 9 .
Figure 9.Comparison between the experimental results and the RVE-3 response for characteristic element length in the range from 2.5 × 10 to 5 × 10 .

Figure 9 .
Figure 9.Comparison between the experimental results and the RVE-3 response for characteristic element length in the range from 2.5 × 10 −4 to 5 × 10 −4 .

Figure 10 .
Figure 10.Comparison between the experimental and the numerical results: (a) transversally loaded unidirectionally reinforced case, and (b) the in-plane shear case.

Figure 10 .
Figure 10.Comparison between the experimental and the numerical results: (a) transversally loaded unidirectionally reinforced case, and (b) the in-plane shear case.

Figure 10 .
Figure 10.Comparison between the experimental and the numerical results: (a) transversally loaded unidirectionally reinforced case, and (b) the in-plane shear case.

Figure 11 .
Figure 11.Damage evolution in the transverse tensile and shear cases: (a) damage index of 0.5 in transverse tension, (b) damage index of 0.75 in transverse tension, (c) failure in transverse tension, (d) damage index of 0.5 in shear, (e) damage index of 0.75 in shear, and (f) failure in shear.

Table 3 .
Dimensions of the additively manufactured specimens.

Table 3 .
Dimensions of the additively manufactured specimens.