The Influence of Filament Orientation on Tensile Stiffness in 3D Printed Structures—Numerical and Experimental Studies

The present study provides a thorough analysis of the influence of filament orientation on the tensile stiffness of 3D-printed structures. This exploration employs a combination of numerical simulations and experimental trials, providing an extensive understanding of additive manufacturing, particularly 3D printing. This process involves layer-by-layer material deposition to produce three-dimensional objects. The examination specifically targets PLA-based 3D printed structures created using Fused Filament Fabrication (FFF) technology and subjects them to rigorous evaluations using a universal tensile testing machine. Additionally, this approach combines Representative Volume Element (RVE) and Classical Lamination Theory (CLT) techniques to extrapolate the mechanical properties of the test material. Although the initial methodology faces challenges in determining the shear modulus with precision, an in-depth investigation results in enhanced accuracy. Furthermore, this study introduces a parametric RVE numerical method, demonstrating its resilience in handling sensitivity to shear modulus. A comparative study of results derived from both the analytical methods and experimental trials involving five series of samples with varied layups reveals that the newly proposed numerical method shows a stronger correlation with the experimental outcomes, delivering a relative error margin of up to 8%.


Introduction
Additive manufacturing (AM), often equated with 3D printing [1] in non-technical contexts [2,3], is a transformative approach to industrial production that employs a layerby-layer addition of materials. The process encompasses a wide array of technologies, including but not limited to binder jetting, material extrusion, and powder bed fusion. The term "additive manufacturing" encapsulates methodologies distinct from traditional subtractive or formative manufacturing techniques. The evolution of AM technologies traces back to the latter half of the 20th century, with significant advancements such as Stereolithography (SLA) [4,5] or Selective Laser Sintering (SLS) [5][6][7]. Specifically, in 1989, Lisa and S. Scott Crump pioneered a material extrusion technique commonly known as Fused Deposition Modeling (FDM) [8,9]. Recognizing its commercial potential, they obtained a patent in 1992 [10] and subsequently founded Stratasys. However, as "Fused Deposition Modeling" and its abbreviation are Stratasys trademarks, the term Fused Filament Fabrication (FFF) has become an alternative descriptor for similar processes [11,12].
The orthotropic nature of FFF printed parts has been emphasized by works like those of Bonada et al. [13], Cuan-Urquizo et al. [14] or Patterson et al. [15,16], who offered multiple approaches of predicting elastic properties of FFF 3D printed parts by investigating the influence of infill pattern on the elastic properties of FFF 3D printed parts or mesoscale modeling. Similar simulations, combined with experimental tests such as those by Fernandez et al., offer predictions on the modulus of elasticity in FFF parts [17]. Additionally,

Manufacturing of the Samples
The initial stage encompassed the printing of the samples. The CAD models of the specimens were initially prepared using SolidWorks 2021 academic edition software, in accordance with [39], and later modified to comply with ASTM D638-14 [37]. The G-code was generated using Ultimaker Cura 5.2.2 software, with key parameters highlighted in Table 1. These parameters were not extensively tested prior to the printing process. Models were produced using an Original Prusa i3 MK3S+ printer equipped with a 0.4 mm brass nozzle. The process entailed the simultaneous manufacture of two specimens positioned flat on the printing bed, as depicted in Figure 1. The selected material for this study was Easy PLA Gray, provided by Fiberlogy (Fiberlab S.A.). The material properties of the material according to the technical data sheet delivered by the producer are given in Table 2 [40]. stresses of experimental measurements of non-standard raster-angle 3D printed samples will be presented and compared with the outcomes of the two numerical methods.

Manufacturing of the Samples
The initial stage encompassed the printing of the samples. The CAD models of the specimens were initially prepared using SolidWorks 2021 academic edition software, in accordance with [39], and later modified to comply with ASTM D638-14 [37]. The G-code was generated using Ultimaker Cura 5.2.2 software, with key parameters highlighted in Table 1. These parameters were not extensively tested prior to the printing process. Models were produced using an Original Prusa i3 MK3S+ printer equipped with a 0.4 mm brass nozzle. The process entailed the simultaneous manufacture of two specimens positioned flat on the printing bed, as depicted in Figure 1. The selected material for this study was Easy PLA Gray, provided by Fiberlogy (Fiberlab S.A.). The material properties of the material according to the technical data sheet delivered by the producer are given in Table  2 [40].

Tensile Testing and Specimen Selection
The printed samples were subsequently categorized and measured before progressing to the second stage of the experimental process. The gauge length section's thickness

Tensile Testing and Specimen Selection
The printed samples were subsequently categorized and measured before progressing to the second stage of the experimental process. The gauge length section's thickness and width of each sample were gauged, as these measurements were required as input for the tensile testing machine.
Tensile testing was performed using a Shimadzu AG-X plus table-top type machine with a maximum load capacity of 50 kN, as depicted in Figure 2a. The initiation of measurements involved supplying the thickness and width values of the samples and setting the number of sets and specimens in each set. An extensometer with a 50-mm base distance was employed to ascertain the difference in the specimen's length before each testing instance. The machine recorded the applied force, including the maximal breaking force, as well as the elongation. A snapshot of a specimen immediately post-testing is illustrated in Figure 2b. and width of each sample were gauged, as these measurements were required as input for the tensile testing machine.
Tensile testing was performed using a Shimadzu AG-X plus table-top type machine with a maximum load capacity of 50 kN, as depicted in Figure 2a. The initiation of measurements involved supplying the thickness and width values of the samples and setting the number of sets and specimens in each set. An extensometer with a 50-mm base distance was employed to ascertain the difference in the specimen's length before each testing instance. The machine recorded the applied force, including the maximal breaking force, as well as the elongation. A snapshot of a specimen immediately post-testing is illustrated in Figure 2b. Before reaching the final series of tensile test measurements, three sets of specimens underwent printing and testing. Several factors, such as specimen dimensions and printing parameters, were evaluated during this preliminary investigation.
Initially, the specimen dimensions were presumed to align with those found in reference [40]. A total of six samples were tested, with two assigned to each filament orientation: 0°, 90°, and 90° × 0°. The latter orientation alternates layers consecutively. The specified degrees correspond to the longitudinal direction of force applied to the specimen during testing or, alternatively, to the longer side of the specimen.
Results from this set revealed that the specimens did not fracture within the gaugelength section, indicating the need for a more in-depth analysis.
The second set of tests encompassed a total of 18 specimens. The dimensions were adjusted to match those of Type 1 and Type 2 test specimens per ASTM standards. For each specimen type, nine were tested, with three dedicated to each filament orientation.
As seen in the attached Figure 3a,b, only specimens with raster angle 90° (specimen numbers 4, 5, 6 and 13, 14, 15) broke in the gauge length section, whereas incorrect breaking remained for all other filament orientations.   Before reaching the final series of tensile test measurements, three sets of specimens underwent printing and testing. Several factors, such as specimen dimensions and printing parameters, were evaluated during this preliminary investigation.
Initially, the specimen dimensions were presumed to align with those found in reference [40]. A total of six samples were tested, with two assigned to each filament orientation: 0 • , 90 • , and 90 • × 0 • . The latter orientation alternates layers consecutively. The specified degrees correspond to the longitudinal direction of force applied to the specimen during testing or, alternatively, to the longer side of the specimen.
Results from this set revealed that the specimens did not fracture within the gaugelength section, indicating the need for a more in-depth analysis.
The second set of tests encompassed a total of 18 specimens. The dimensions were adjusted to match those of Type 1 and Type 2 test specimens per ASTM standards. For each specimen type, nine were tested, with three dedicated to each filament orientation.
As seen in the attached Figure 3a,b, only specimens with raster angle 90 • (specimen numbers 4, 5, 6 and 13, 14, 15) broke in the gauge length section, whereas incorrect breaking remained for all other filament orientations.
and width of each sample were gauged, as these measurements were required as input for the tensile testing machine.
Tensile testing was performed using a Shimadzu AG-X plus table-top type machine with a maximum load capacity of 50 kN, as depicted in Figure 2a. The initiation of measurements involved supplying the thickness and width values of the samples and setting the number of sets and specimens in each set. An extensometer with a 50-mm base distance was employed to ascertain the difference in the specimen's length before each testing instance. The machine recorded the applied force, including the maximal breaking force, as well as the elongation. A snapshot of a specimen immediately post-testing is illustrated in Figure 2b. Before reaching the final series of tensile test measurements, three sets of specimens underwent printing and testing. Several factors, such as specimen dimensions and printing parameters, were evaluated during this preliminary investigation.
Initially, the specimen dimensions were presumed to align with those found in reference [40]. A total of six samples were tested, with two assigned to each filament orientation: 0°, 90°, and 90° × 0°. The latter orientation alternates layers consecutively. The specified degrees correspond to the longitudinal direction of force applied to the specimen during testing or, alternatively, to the longer side of the specimen.
Results from this set revealed that the specimens did not fracture within the gaugelength section, indicating the need for a more in-depth analysis.
The second set of tests encompassed a total of 18 specimens. The dimensions were adjusted to match those of Type 1 and Type 2 test specimens per ASTM standards. For each specimen type, nine were tested, with three dedicated to each filament orientation.
As seen in the attached Figure 3a,b, only specimens with raster angle 90° (specimen numbers 4, 5, 6 and 13, 14, 15) broke in the gauge length section, whereas incorrect breaking remained for all other filament orientations.   This realization spurred another iteration in the specimen preparation process. For Type 1 specimens, the CAD model was adjusted so that the curves of the specimen's longer edges were modeled as a single entity, employing the Spline function rather than connecting two arcs and a straight line, as depicted in Figure 4a. Additionally, a new specimen type was tested, which featured "V" shapes in place of the arcs connecting the gauge-length section of the specimen with its wider section, as shown in Figure 4b. Furthermore, during the G-code preparation, the layer width was fine-tuned from 0.4 mm to 0.2 mm. Consequently, the third set consisted of 30 specimens: 18 "V"-shaped specimens with three filament orientations each, varying in layer width (0.2 mm and 0.4 mm), with three per setting; and 12 specimens implementing the Spline function in the CAD model, organized similarly to the "V"-shaped specimens, albeit excluding the 90 • orientation since it had already demonstrated proper breakage in the second set of measurements.
This realization spurred another iteration in the specimen preparation process. For Type 1 specimens, the CAD model was adjusted so that the curves of the specimen's longer edges were modeled as a single entity, employing the Spline function rather than connecting two arcs and a straight line, as depicted in Figure 4a. Additionally, a new specimen type was tested, which featured "V" shapes in place of the arcs connecting the gauge-length section of the specimen with its wider section, as shown in Figure 4b. Furthermore, during the G-code preparation, the layer width was fine-tuned from 0.4 mm to 0.2 mm. Consequently, the third set consisted of 30 specimens: 18 "V"-shaped specimens with three filament orientations each, varying in layer width (0.2 mm and 0.4 mm), with three per setting; and 12 specimens implementing the Spline function in the CAD model, organized similarly to the "V"-shaped specimens, albeit excluding the 90° orientation since it had already demonstrated proper breakage in the second set of measurements.  Despite these alterations, the set of specimens did not yield any significant changes in the breaking location. Most specimens fractured either outside or at the end of the gauge-length section. Figure 5 displays specimens from the third measurement set that broke at the correct location. Their indices are specified at the top of the specimens, with layer width and filament orientation (090 refers to 90° × 0°) notated at the bottom. Therefore, Type 2 specimens in Spline variation have been selected for further testing.  Despite these alterations, the set of specimens did not yield any significant changes in the breaking location. Most specimens fractured either outside or at the end of the gaugelength section. Figure 5 displays specimens from the third measurement set that broke at the correct location. Their indices are specified at the top of the specimens, with layer width and filament orientation (090 refers to 90 • × 0 • ) notated at the bottom. Therefore, Type 2 specimens in Spline variation have been selected for further testing.
This realization spurred another iteration in the specimen preparation process. For Type 1 specimens, the CAD model was adjusted so that the curves of the specimen's longer edges were modeled as a single entity, employing the Spline function rather than connecting two arcs and a straight line, as depicted in Figure 4a. Additionally, a new specimen type was tested, which featured "V" shapes in place of the arcs connecting the gauge-length section of the specimen with its wider section, as shown in Figure 4b. Furthermore, during the G-code preparation, the layer width was fine-tuned from 0.4 mm to 0.2 mm. Consequently, the third set consisted of 30 specimens: 18 "V"-shaped specimens with three filament orientations each, varying in layer width (0.2 mm and 0.4 mm), with three per setting; and 12 specimens implementing the Spline function in the CAD model, organized similarly to the "V"-shaped specimens, albeit excluding the 90° orientation since it had already demonstrated proper breakage in the second set of measurements.  Despite these alterations, the set of specimens did not yield any significant changes in the breaking location. Most specimens fractured either outside or at the end of the gauge-length section. Figure 5 displays specimens from the third measurement set that broke at the correct location. Their indices are specified at the top of the specimens, with layer width and filament orientation (090 refers to 90° × 0°) notated at the bottom. Therefore, Type 2 specimens in Spline variation have been selected for further testing.  A stress-strain curve was prepared for each sample based on the extensometer results to determine Young's modulus from sample measurements. Young's modulus was derived directly from the slope in the linear elastic region. The highest stress value represented the failure load, or the stress at breakage. For instance, Figure 6 illustrates a typical stress-strain curve for a sample with a 0 • raster angle. The black line symbolizes Young's modulus slope, equating to 3600 MPa, calculated as the difference in stress values of 18 and 3.6 and strain values of 0.005 and 0.001. The pink dashed line portrays the failure load, determined from the highest point of the stress-strain curve and transferred to the graph's vertical axis, resulting in a value of 56 MPa. Similar calculations were made for all other specimens, following the guidelines outlined in the standard [37]. Table 3 summarizes the results from all three sets of preliminary studies.
A stress-strain curve was prepared for each sample based on the extensometer results to determine Young's modulus from sample measurements. Young's modulus was derived directly from the slope in the linear elastic region. The highest stress value represented the failure load, or the stress at breakage. For instance, Figure 6 illustrates a typical stress-strain curve for a sample with a 0° raster angle. The black line symbolizes Young's modulus slope, equating to 3600 MPa, calculated as the difference in stress values of 18 and 3.6 and strain values of 0.005 and 0.001. The pink dashed line portrays the failure load, determined from the highest point of the stress-strain curve and transferred to the graph's vertical axis, resulting in a value of 56 MPa. Similar calculations were made for all other specimens, following the guidelines outlined in the standard [37]. Table 3 summarizes the results from all three sets of preliminary studies.  Calculations were made for Young's modulus and failure stress for all specimens in each set. Although relatively high standard deviations of the values have been achieved due to various types of tested samples, a strong correlation was observed between these two properties. The values obtained from the calculations align well with the technical data sheet for the filament [16]. The reference provides Young's modulus and tensile strength at break of 3500 MPa and 53 MPa, respectively, which closely match the values for a 0° filament orientation. This raster angle is closest in strength to isotropic PLA.

Microscopic Testing
Geometry preparation encompassed an examination of the specimen's microstructure post-breakage. To aid in this investigation, a specimen with a 0° filament orientation was inspected under a microscope.
The microstructure presented in Figure 7 reveals that the gaps between filament strips are not uniformly homogeneous. The magnitude of these gaps influences the overall strength of a 3D-printed material. Consequently, the largest gap was selected and used as an input for the Representative Volume Element (RVE) model. Figure 7 displays a red square with dashed lines that approximates the cross-sectional area of the RVE model. The overall dimensions of this area align with the layer size, specifically 0.2 mm in height and 0.4 mm in width, in accordance with Figure 7. The exact geometry of the RVE will be discussed more thoroughly in the following sections of this article.  Calculations were made for Young's modulus and failure stress for all specimens in each set. Although relatively high standard deviations of the values have been achieved due to various types of tested samples, a strong correlation was observed between these two properties. The values obtained from the calculations align well with the technical data sheet for the filament [16]. The reference provides Young's modulus and tensile strength at break of 3500 MPa and 53 MPa, respectively, which closely match the values for a 0 • filament orientation. This raster angle is closest in strength to isotropic PLA.

Microscopic Testing
Geometry preparation encompassed an examination of the specimen's microstructure post-breakage. To aid in this investigation, a specimen with a 0 • filament orientation was inspected under a microscope.
The microstructure presented in Figure 7 reveals that the gaps between filament strips are not uniformly homogeneous. The magnitude of these gaps influences the overall strength of a 3D-printed material. Consequently, the largest gap was selected and used as an input for the Representative Volume Element (RVE) model. Figure 7 displays a red square with dashed lines that approximates the cross-sectional area of the RVE model. The overall dimensions of this area align with the layer size, specifically 0.2 mm in height and 0.4 mm in width, in accordance with Figure 7. The exact geometry of the RVE will be discussed more thoroughly in the following sections of this article.

Numerical and Analytical Studies
This section presents two approaches to predicting the elastic properties (Young's moduli) of analyzed 3D printed structures with different raster angle layups. The first method combines Representative Volume Element (RVE) and composite materials using Classical Lamination Theory (CLT). In the second approach, RVE modeling is performed across the entire spectrum of raster angles, and the elastic properties of multi-layered 3D prints are found as the weighted average of the properties of individual layers inclined at different angles.

Mixed RVE and CLT Approach
In this approach, based on the microscopic image of the structure, a finite element model of the unit cell is created, and material properties in the main orthotropic directions are determined. Later, composite material CLT is applied to find the elastic properties of multi-layered 3D prints with different raster angle arrangements.

Representative Volume Element Approach
In the initial step, a finite element model of a unit volume of the structure is created based on the microscopic image presented in Figure 7. The model was prepared in Space-Claim, ANSYS 2022 R1 in-built modeling software, and is presented in Figure 8a,b. As an input material, isotropic PLA properties were used, provided by the manufacturer of the filament [41], specifically: Young modulus E = 3500 MPa, whereas Poisson Ratio = 0.35, based on the literature [42,43]. Later, the model was discretized using SOLID186 hexahedral elements [44], with an element size of 0.019 mm, a total number of elements equal to 6615, and 30,681 nodes-see Figure 8c.

Numerical and Analytical Studies
This section presents two approaches to predicting the elastic properties (Young's moduli) of analyzed 3D printed structures with different raster angle layups. The first method combines Representative Volume Element (RVE) and composite materials using Classical Lamination Theory (CLT). In the second approach, RVE modeling is performed across the entire spectrum of raster angles, and the elastic properties of multi-layered 3D prints are found as the weighted average of the properties of individual layers inclined at different angles.

Mixed RVE and CLT Approach
In this approach, based on the microscopic image of the structure, a finite element model of the unit cell is created, and material properties in the main orthotropic directions are determined. Later, composite material CLT is applied to find the elastic properties of multi-layered 3D prints with different raster angle arrangements.

Representative Volume Element Approach
In the initial step, a finite element model of a unit volume of the structure is created based on the microscopic image presented in Figure 7. The model was prepared in Space-Claim, ANSYS 2022 R1 in-built modeling software, and is presented in Figure 8a,b. As an input material, isotropic PLA properties were used, provided by the manufacturer of the filament [41], specifically: Young modulus E = 3500 MPa, whereas Poisson Ratio ν = 0.35, based on the literature [42,43]. Later, the model was discretized using SOLID186 hexahedral elements [44], with an element size of 0.019 mm, a total number of elements equal to 6615, and 30,681 nodes-see Figure 8c.

Numerical and Analytical Studies
This section presents two approaches to predicting the elastic properties (Young's moduli) of analyzed 3D printed structures with different raster angle layups. The first method combines Representative Volume Element (RVE) and composite materials using Classical Lamination Theory (CLT). In the second approach, RVE modeling is performed across the entire spectrum of raster angles, and the elastic properties of multi-layered 3D prints are found as the weighted average of the properties of individual layers inclined at different angles.

Mixed RVE and CLT Approach
In this approach, based on the microscopic image of the structure, a finite element model of the unit cell is created, and material properties in the main orthotropic directions are determined. Later, composite material CLT is applied to find the elastic properties of multi-layered 3D prints with different raster angle arrangements.

Representative Volume Element Approach
In the initial step, a finite element model of a unit volume of the structure is created based on the microscopic image presented in Figure 7. The model was prepared in Space-Claim, ANSYS 2022 R1 in-built modeling software, and is presented in Figure 8a,b. As an input material, isotropic PLA properties were used, provided by the manufacturer of the filament [41], specifically: Young modulus E = 3500 MPa, whereas Poisson Ratio = 0.35, based on the literature [42,43]. Later, the model was discretized using SOLID186 hexahedral elements [44], with an element size of 0.019 mm, a total number of elements equal to 6615, and 30,681 nodes-see Figure 8c.  In order to predict all necessary mechanical properties, six simulations with varying load types and directions are performed. The first three simulations were used to determine Young's moduli in three directions: E x , E y , E z , and Poisson's ratios ν yz , ν xz , and ν xysee Figure 9a-c, respectively. The remaining simulations were used to determine shear (Kirchhoff's) moduli G yz , G xz , and G xy -see Figure 9d-f, respectively. A schematic of the simulation boundary conditions is presented in Figure 10a, supplied with views from Ansys 2022 R1 software of exemplary constraints in tensile and shear simulations in Figure 10b-d. In order to predict all necessary mechanical properties, six simulations with varying load types and directions are performed. The first three simulations were used to determine Young's moduli in three directions: Ex, Ey, Ez, and Poisson's ratios νyz, νxz, and νxysee Figure 9a-c, respectively. The remaining simulations were used to determine shear (Kirchhoff's) moduli Gyz, Gxz, and Gxy-see Figure 9d-f, respectively. A schematic of the simulation boundary conditions is presented in Figure 10a, supplied with views from Ansys 2022 R1 software of exemplary constraints in tensile and shear simulations in Figure  10b-d.   In order to predict all necessary mechanical properties, six simulations with varying load types and directions are performed. The first three simulations were used to determine Young's moduli in three directions: Ex, Ey, Ez, and Poisson's ratios νyz, νxz, and νxysee Figure 9a-c, respectively. The remaining simulations were used to determine shear (Kirchhoff's) moduli Gyz, Gxz, and Gxy-see Figure 9d-f, respectively. A schematic of the simulation boundary conditions is presented in Figure 10a, supplied with views from Ansys 2022 R1 software of exemplary constraints in tensile and shear simulations in Figure  10b-d.   Tensile simulations were conducted according to the schematic shown in Figure 10a. Three displacement boundary conditions were applied to the RVE model: a surface parallel to the one with force applied; a line and a node, each constraining one degree of freedom, as shown in Figure 10b. When it comes to shear simulations, additional APDL code is needed in order to prevent undesirable deformations. In all three cases, the surface with applied force was coupled in two directions: one parallel and one perpendicular to the direction of the force vector. As shown in Figure 10c,d, force is applied to the surface "front", whereas u z and u x are directions of node coupling of the surface.
Having prepared and carried out the simulations, results were gathered and used to calculate all nine mechanical properties. Exemplary deformation plots extracted from Ansys software are shown in Figure 11, whereas calculation results are shown in Table 4. Tensile simulations were conducted according to the schematic shown in Figure 10a. Three displacement boundary conditions were applied to the RVE model: a surface parallel to the one with force applied; a line and a node, each constraining one degree of freedom, as shown in Figure 10b. When it comes to shear simulations, additional APDL code is needed in order to prevent undesirable deformations. In all three cases, the surface with applied force was coupled in two directions: one parallel and one perpendicular to the direction of the force vector. As shown in Figure 10c,d, force is applied to the surface "front", whereas uz and ux are directions of node coupling of the surface.
Having prepared and carried out the simulations, results were gathered and used to calculate all nine mechanical properties. Exemplary deformation plots extracted from Ansys software are shown in Figure 11, whereas calculation results are shown in Table 4.  The above results were used as input into an analytical method explained below.

Composites Classical Lamination Theory-Prediction of Young's Modulus
According to the theory described in [45], Classical Lamination Theory (CLT) is an approach allowing the determination of various mechanical properties of an orthotropic composite material. A laminate is an organized pile of uni-directional, bonded plies (or layers) of composite materials with a given direction of its fibers. In addition, the layers have much bigger lengths than thicknesses. Since orthotropic material is considered, the direction of stiffness, i.e., Young's moduli, needs to be taken into account: where 1 and 2 subscripts refer to the direction of a property in relation to the direction of the plies, parallel and perpendicular, respectively. Analogically, shear forces are defined using shear (Kirchhoff's) modulus, which links shear stress with shear strain :  The above results were used as input into an analytical method explained below.

Composites Classical Lamination Theory-Prediction of Young's Modulus
According to the theory described in [45], Classical Lamination Theory (CLT) is an approach allowing the determination of various mechanical properties of an orthotropic composite material. A laminate is an organized pile of uni-directional, bonded plies (or layers) of composite materials with a given direction of its fibers. In addition, the layers have much bigger lengths than thicknesses. Since orthotropic material is considered, the direction of stiffness, i.e., Young's moduli, needs to be taken into account: where 1 and 2 subscripts refer to the direction of a property in relation to the direction of the plies, parallel and perpendicular, respectively. Analogically, shear forces are defined using shear (Kirchhoff's) modulus, which links shear stress τ with shear strain γ: In the case of plane stress, namely, a plate undergoing stresses in various directions within a plane, Poisson's ratios need to be considered: Using all the properties mentioned above, a matrix equation can be derived: The 3 × 3 matrix, called the compliance matrix, consists of various Q's, which are reduced stiffnesses: To transform the stresses and strains into the principal material directions, the matrix takes the form: where x and y are axes corresponding to the direction of the applied load. Components in the lamina stress matrix for a ply with fibers oriented at θ may be represented as: Q 22 = Q 11 n 4 + 2(Q 12 + 2Q 66 )m 2 n 2 + Q 22 m 4 (12) where m = cos θ and n = sin θ. It can be seen that the resultant stress will produce both normal and shear stress on the fibers of the laminate, even in the case of a normal load applied to the whole laminate.
Matrices denoted as A, B, and D are introduced to arrange all possible loads acting on the laminate. They correspond to the thickness of the ply, its distance from the midplane of the laminate, and specific loads. Hence, the A matrix, namely, the extensional stiffness matrix, relates to normal stresses and strains. Taking into account the geometry of the layers and the whole laminate in this research, the second term of the matrix shall be represented as t k : Similarly, the coupling stiffness matrix, namely the B matrix, relates bending strains to normal stress, or vice versa: Eventually, the D matrix, which is called the bending stiffness matrix, corresponds to the number of plate curvatures and the bending moments: Young's Modulus E x , can be determined using components from all three matrices mentioned above and strain in the x direction: Young's Modulus in the y direction E y , and shear modulus G xy , can be found similarly with the use of terms from the A, B, and D matrices: In order to derive an equation for Poisson's Ratios, strain in the y direction while stress is applied in the x direction is needed: Then, ν xy will be expressed as a ratio of strains in the x and y directions, respectively: The equations presented above were used to determine all necessary mechanical properties using data obtained from the numerical simulations of the RVE model.

Angular RVE-Raster Angle Parametrized Model
A second numerical approach was also proposed to supplement the methods mentioned earlier. Based on the geometry presented in Section 3.1.1, a parameterized simulation was prepared to determine the behavior of the RVE's mechanical properties across the entire spectrum of raster angles.
The geometry remained unchanged; nevertheless, the mesh was prepared again. A SOLID187 tetrahedral element with midside nodes was used in this case [44], achieving 90,194 elements and 134,355 nodes. The mesh is depicted in Figure 12  The parameter in the simulation was set directly in the model's geometry. It is based on rotation around a vertical axis, as seen in Figure 12. The step size was set to 2.5°, which returned 36 × 6 = 216 simulations. During the simulations, necessary data was preserved, namely: three surface areas (relating to the front, top, and side walls of the model) and displacements; in tensile simulations, one per each axis, nine in total; in shear simulations, one per each simulation, three in total.
All this data was necessary for calculating the averaged Young's modulus, carried out in the second step. In order to do that, one has to familiarize with the methodology explained below.
The specimens tested in this study consist of 16 layers of nearly homogenous filament strips. Initially, one filament layer will be considered. From Hooke's Law: where F is the tensile force acting on the layer; l is the initial length of the layer; w is the width of the layer; t is the thickness of the layer; and ∆ is the displacement of the layer Solving for the force in Equation (27) produces: One can translate this equation to several layers, representing it as a sum of forces per layer, even though the force magnitude at a specific layer is unknown.   The parameter in the simulation was set directly in the model's geometry. It is based on rotation around a vertical axis, as seen in Figure 12. The step size was set to 2.5°, which returned 36 × 6 = 216 simulations. During the simulations, necessary data was preserved, namely: three surface areas (relating to the front, top, and side walls of the model) and displacements; in tensile simulations, one per each axis, nine in total; in shear simulations, one per each simulation, three in total.
All this data was necessary for calculating the averaged Young's modulus, carried out in the second step. In order to do that, one has to familiarize with the methodology explained below.
The specimens tested in this study consist of 16 layers of nearly homogenous filament strips. Initially, one filament layer will be considered. From Hooke's Law: where F is the tensile force acting on the layer; l is the initial length of the layer; w is the width of the layer; t is the thickness of the layer; and ∆ is the displacement of the layer Solving for the force in Equation (27) produces: One can translate this equation to several layers, representing it as a sum of forces per layer, even though the force magnitude at a specific layer is unknown. The parameter in the simulation was set directly in the model's geometry. It is based on rotation around a vertical axis, as seen in Figure 12. The step size was set to 2.5 • , which returned 36 × 6 = 216 simulations. During the simulations, necessary data was preserved, namely: three surface areas (relating to the front, top, and side walls of the model) and displacements; in tensile simulations, one per each axis, nine in total; in shear simulations, one per each simulation, three in total.
All this data was necessary for calculating the averaged Young's modulus, carried out in the second step. In order to do that, one has to familiarize with the methodology explained below.
The specimens tested in this study consist of 16 layers of nearly homogenous filament strips. Initially, one filament layer will be considered. From Hooke's Law: where F is the tensile force acting on the layer; l is the initial length of the layer; w is the width of the layer; t is the thickness of the layer; and ∆l is the displacement of the layer. Solving for the force in Equation (27) produces: One can translate this equation to several layers, representing it as a sum of forces per layer, even though the force magnitude at a specific layer is unknown.
where n is the number of layers. Implementing Equation (28) into Equation (29), while bearing in mind that all components, in addition to Young's modulus and the layer thickness, of Equation (28) are constant for all layers, Equation (29) can be represented as a sum of Young's moduli multiplied by consecutive thicknesses: The technique is supplied with a schematic shown in Figure 14. Using such an approach, both orthotropic Young's moduli E 1 and E 2 , were calculated based on the data obtained from simulations of angular RVE.  (29) where n is the number of layers. Implementing Equation (28) into Equation (29), while bearing in mind that all components, in addition to Young's modulus and the layer thickness, of Equation (28) are constant for all layers, Equation (29) can be represented as a sum of Young's moduli multiplied by consecutive thicknesses: The technique is supplied with a schematic shown in Figure 14. Using such an approach, both orthotropic Young's moduli 1 and 2 , were calculated based on the data obtained from simulations of angular RVE.

Comparison of the Models
In this section, calculations results for the whole spectrum of raster angles of the two methods, RVE-CLT and angular RVE, will be compared based on plots from Figure 15.

Comparison of the Models
In this section, calculations results for the whole spectrum of raster angles of the two methods, RVE-CLT and angular RVE, will be compared based on plots from Figure 15.  (29) where n is the number of layers. Implementing Equation (28) into Equation (29), while bearing in mind that all components, in addition to Young's modulus and the layer thickness, of Equation (28) are constant for all layers, Equation (29) can be represented as a sum of Young's moduli multiplied by consecutive thicknesses: The technique is supplied with a schematic shown in Figure 14. Using such an approach, both orthotropic Young's moduli 1 and 2 , were calculated based on the data obtained from simulations of angular RVE.

Comparison of the Models
In this section, calculations results for the whole spectrum of raster angles of the two methods, RVE-CLT and angular RVE, will be compared based on plots from Figure 15.  The above graphs show the results of Young's moduli, shear modulus, and Poisson ratio calculated for the two methods presented earlier in the study. Both methods, when analyzed separately, preserve a specific behavior; for instance, in the case of RVE-CLT, for both Young moduli E 1 and E 2 , lowest raster angles are for 45 • and −45 • , whereas in the vicinity of the same raster angles, plots of angular RVE depict either an increase or decrease. Nevertheless, more detailed conclusions cannot be drawn from either of the four pairs of plots, as they do not depict any similar behavior.

RVE-CLT Approach-Troubleshooting
As shown in Figure 15, it is visible that the mixed RVE-CLT approach significantly deviates from values of Young's moduli in the case of samples with layups different than [0,0] and [90,90]. Troubleshooting of this issue will be performed based on the workflow analysis of this approach. First, it is worth noticing that the input parameters for the RVE finite element model are isotropic, elastic properties of the material (Young's modulus E and Poisson's ratio ν). The output of RVE modeling used as an input to Classical Lamination Theory are Young's moduli in two orthotropic directions, E 1 and E 2 , shear modulus G 12 , and Poisson's ratio ν 12 . Bearing in mind the results of the preliminary study, it is worth noting that the result of E 2 , obtained from RVE (E 2,RVE = 2938 MPa), can be treated as expected. It means that potential failure in the approach can be hidden in the determination of shear modulus G 12 or Poisson's ratio ν 12 . On the other hand, Poisson's ratio is a relatively simple quantity to determine as a fraction of lateral shortening to elongation of the sample. Hence, it suggests that a solution should be sought to determine the shear modulus.
Equation (31) supplied with a schematic from Figure 16 defines the shear modulus as a ratio of shear stress and strain. On the other hand, for isotropic materials, the relationship between Young's modulus, Poisson's ratio, and shear modulus is given as [46]: The above graphs show the results of Young's moduli, shear modulus, and Poisson ratio calculated for the two methods presented earlier in the study. Both methods, when analyzed separately, preserve a specific behavior; for instance, in the case of RVE-CLT, for both Young moduli 1 and 2 , lowest raster angles are for 45° and −45°, whereas in the vicinity of the same raster angles, plots of angular RVE depict either an increase or decrease. Nevertheless, more detailed conclusions cannot be drawn from either of the four pairs of plots, as they do not depict any similar behavior.

RVE-CLT Approach-Troubleshooting
As shown in Figure 15, it is visible that the mixed RVE-CLT approach significantly deviates from values of Young's moduli in the case of samples with layups different than [0,0] and [90,90]. Troubleshooting of this issue will be performed based on the workflow analysis of this approach. First, it is worth noticing that the input parameters for the RVE finite element model are isotropic, elastic properties of the material (Young's modulus E and Poisson's ratio ν). The output of RVE modeling used as an input to Classical Lamination Theory are Young's moduli in two orthotropic directions, E1 and E2, shear modulus G12, and Poisson's ratio ν12. Bearing in mind the results of the preliminary study, it is worth noting that the result of E2, obtained from RVE (E2,RVE = 2938 MPa), can be treated as expected. It means that potential failure in the approach can be hidden in the determination of shear modulus G12 or Poisson's ratio ν12. On the other hand, Poisson's ratio is a relatively simple quantity to determine as a fraction of lateral shortening to elongation of the sample. Hence, it suggests that a solution should be sought to determine the shear modulus.
Equation (31) supplied with a schematic from Figure 16 defines the shear modulus as a ratio of shear stress and strain. On the other hand, for isotropic materials, the relationship between Young's modulus, Poisson's ratio, and shear modulus is given as [46]: To see if the numerical model is constructed correctly, the RVE model will be prepared for a well-known material-structural steel. A FE model of shearing (shear force F = 100 N) of the cube of 10 mm × 10 mm × 10 mm made from steel with E = 200 GPa and ν = 0.3 is prepared and presented in Figure 17a,b. To see if the numerical model is constructed correctly, the RVE model will be prepared for a well-known material-structural steel. A FE model of shearing (shear force F = 100 N) of the cube of 10 mm × 10 mm × 10 mm made from steel with E = 200 GPa and ν = 0.3 is prepared and presented in Figure 17a On the other hand, according to Equation (32), the shear modulus for steel is equal to: This discrepancy indicates that the numerical model must be carefully reviewed. To do it, in-plane shear stress is displayed: Figure 18 shows that the shear stress in the middle of the wall is different than on its boundaries, which can be explained by Saint-Venant's principle [46]. Taking the value of shear stress from the middle of the wall and substituting it with the formula for shear modulus:  Based on this model, the shear modulus is equal to: On the other hand, according to Equation (32), the shear modulus for steel is equal to: This discrepancy indicates that the numerical model must be carefully reviewed. To do it, in-plane shear stress is displayed: Figure 18 shows that the shear stress in the middle of the wall is different than on its boundaries, which can be explained by Saint-Venant's principle [46]. Taking the value of shear stress from the middle of the wall and substituting it with the formula for shear modulus: is an acceptable result. Applying this observation to the 3D printed structures analyzed in this study, one obtains: Since the RVE model is assumed to be isotropic, it is possible to use the analytical formula: Substituting these results to the Classical Lamination Theory and confronting them with the angular RVE approach, one obtains updated graphs of orthotropic material properties.
Similarly to Figure 15, the graphs shown in Figure 19 represent the results of Young's moduli, shear modulus, and Poisson ratio calculated, however, with an altered value of shear modulus. In this case, the RVE angular and RVE-CLT plots correspond much better. Full-spectrum plots of Young's moduli maintain similar behavior, with either increase or decrease in the vicinity of 45 • or −45 • raster angles, whereas interchanging peaks at 90 • , −90 • and 0, for instance E 1 attains its maximum value at 0 • , on the other hand, at 0 • , E 2 attains its minimum. Shear modulus and Poisson's ratio plots correspond slightly worse, yet after changing the methodology for shear modulus determination, the situation is improved, especially in the case of Poisson's ratio.
On the other hand, according to Equation (32), the shear modulus for steel is equal to: This discrepancy indicates that the numerical model must be carefully reviewed. To do it, in-plane shear stress is displayed: Figure 18 shows that the shear stress in the middle of the wall is different than on its boundaries, which can be explained by Saint-Venant's principle [46]. Taking the value of shear stress from the middle of the wall and substituting it with the formula for shear modulus: Figure 18. In-plane shear resultant stress of a steel cube.
is an acceptable result. Applying this observation to the 3D printed structures analyzed in this study, one obtains: Since the RVE model is assumed to be isotropic, it is possible to use the analytical formula: Substituting these results to the Classical Lamination Theory and confronting them with the angular RVE approach, one obtains updated graphs of orthotropic material properties.
Similarly to Figure 15, the graphs shown in Figure 19 represent the results of Young's moduli, shear modulus, and Poisson ratio calculated, however, with an altered value of shear modulus. In this case, the RVE angular and RVE-CLT plots correspond much better. Full-spectrum plots of Young's moduli maintain similar behavior, with either increase or decrease in the vicinity of 45° or −45° raster angles, whereas interchanging peaks at 90°, −90° and 0, for instance 1 attains its maximum value at 0°, on the other hand, at 0°, 2 attains its minimum. Shear modulus and Poisson's ratio plots correspond slightly worse, yet after changing the methodology for shear modulus determination, the situation is improved, especially in the case of Poisson's ratio.

Results
In this section, the results of experimental studies will be presented, including stressstrain curves, averaged Young's moduli, and failure stresses. In the final stage, stiffnesses will be confronted with predictions obtained from an adjusted RVE-CLT methodology as well as the RVE method with a raster angle parametrized model (angular RVE).

Experimental Results
Bearing in mind the preliminary tests, a set containing 50 specimens was 3D printed, preserving all machine settings, dimensions, CAD models, and printing process settings. Since the layer height was set to 0.2 mm and the specimens were 3.2 mm thick, they consisted of 16 layers. Five series, 10 specimens each, of different layups were examined: It means that in the case of the first two layups, the layers maintained the same raster angle for all 16 layers. Series S3, S4, and S5 included interchanging layers with a symmetry plane between the 8th and 9th layers, bearing in mind that the raster angle is measured according to the symmetry axis of the specimen, or, in other words, the longer side of the sample. Conversely, other layups displayed significantly higher strains at failure, ranging from 0.03 to 0.04, compared to a maximum value of 0.02 in S1. This could be attributed to the interchanging layers of these samples. One plot that significantly deviates from the others in terms of strain is the [−45/45] layup. For most samples, the strain exceeded 0.04, which is markedly higher than other layups. Normal and shear forces in this particular arrangement led to the largest plastic deformation. A notable distinction between S1 and S2, compared to the other series, is the uniformity in the samples' behavior under tensile stress. Plots of individual results, indicated by thin lines on the graphs, demonstrate that failure occurred over a broader range of strains in the S3, S4, and S5 sets compared to the first two.
When comparing all average results of stress and strain relationships, as displayed in Figure 22, a nearly pure brittle fracture was observed in the case of 0 • filament orientation. However, increased brittleness was impacted by the de-bonding of the filament strips in a sample with a [90/90] layup, as the force direction was perpendicular to the orientation of the filament strips.
Average Young's moduli and failure stresses are documented in Figures 21 and 23. The highest Young's modulus was obtained for S1, equating to 3542 ± 51 MPa, with a decrease observed for increasing filament angles, reaching 2726 ± 81 MPa in the [90,90] layup. This also correlates to the stress vs. strain plots. The distribution of Young's moduli magnitudes relates to the highest measured stresses. Although the angular difference between the 0 • and 30 • filament orientations is the same as for the last two pairs, the difference in Young's moduli is much higher for the first pair, approximately 300 MPa compared to 150 MPa. This could be due to the type of stress experienced in the [0,0] filament orientation. A specimen subjected to pure axial stress, where shear stress is absent, results in a significantly higher strength than any other orientation. In addition to pure tensile (0 • ) and shear (90 • ) loads, the laminate is subjected to both. Given that interlayer bonding in the 90 • orientation needs to be taken into account, the overall strength will be considerably lower.
which is markedly higher than other layups. Normal and shear forces in this particular arrangement led to the largest plastic deformation. A notable distinction between S1 and S2, compared to the other series, is the uniformity in the samples' behavior under tensile stress. Plots of individual results, indicated by thin lines on the graphs, demonstrate that failure occurred over a broader range of strains in the S3, S4, and S5 sets compared to the first two.     When comparing all average results of stress and strain relationships, as displayed in Figure 22, a nearly pure brittle fracture was observed in the case of 0° filament orientation. However, increased brittleness was impacted by the de-bonding of the filament    When comparing all average results of stress and strain relationships, as displayed in Figure 22, a nearly pure brittle fracture was observed in the case of 0° filament orientation. However, increased brittleness was impacted by the de-bonding of the filament In all instances, failure stress corresponds to the obtained Young's modulus, with insignificant or very low standard deviations. Again, as per the graphs, specimens with a 0 • orientation failed under the highest stress, amounting to 57 ± 1 MPa, while 90 • specimens failed the quickest-under only 27 ± 1 MPa. A range of fracture types were observed in the study.

Experimental vs. Numerical Results-Comparison
As the final stage in this study, a graph with plots of averaged Young's moduli with all methodologies included is shown in Figure 24.
The strength behavior along the spectrum of raster angles, namely lower raster angles, as in S4 or S5, differs less in strength than higher raster angles, as in S1 or S2. In the S3 and S4 series, experimental results correlate well with angular RVE results: 3000 MPa to 3100 MPa and 2900 MPa to 3000 MPa-in both cases, relative error maintained a lower level than 3.5%. Nevertheless, the least indefinite results between specific methods can be spotted for the S1 and S2 series, with the least approximated results for a 0 • filament orientation: 3500 (experimental) and 3500 (RVE-CLT and angular RVE).
0° orientation failed under the highest stress, amounting to 57 ± 1 MPa, while 90° speci-mens failed the quickest-under only 27 ± 1 MPa. A range of fracture types were observed in the study.

Experimental vs. Numerical Results-Comparison
As the final stage in this study, a graph with plots of averaged Young's moduli with all methodologies included is shown in Figure 24. The strength behavior along the spectrum of raster angles, namely lower raster angles, as in S4 or S5, differs less in strength than higher raster angles, as in S1 or S2. In the S3 and S4 series, experimental results correlate well with angular RVE results: 3000 MPa to 3100 MPa and 2900 MPa to 3000 MPa-in both cases, relative error maintained a lower level than 3.5%. Nevertheless, the least indefinite results between specific methods can be spotted for the S1 and S2 series, with the least approximated results for a 0° filament orientation: 3500 (experimental) and 3500 (RVE-CLT and angular RVE).
An angular RVE plot undoubtedly reflects experimental results better than an RVE-CLT plot. CLT is a methodology developed with the intention of utilizing it with composite materials. The most notable difference is layer bonding. The 3D-printed samples consist of gaps, whereas composite laminates do not. Hence, this discrepancy in results may An angular RVE plot undoubtedly reflects experimental results better than an RVE-CLT plot. CLT is a methodology developed with the intention of utilizing it with composite materials. The most notable difference is layer bonding. The 3D-printed samples consist of gaps, whereas composite laminates do not. Hence, this discrepancy in results may come from the fact that FFF 3D printed material does not strictly comply with assumptions made for Classical Lamination Theory regarding composite materials.

Discussion
Before reaching numerical procedures, 3D-printed specimens were manufactured and tested using a tensile test machine. Preliminary tests proved to be necessary before conducting experiments. They helped with attaining a suitable way to produce specimens that could be used during the measurements. Applying dimensions from ASTM standards allowed for proper breaking in the case of [90,90] layup. In addition, calculations based on data obtained from measurements complied with assumptions made before attempting the study, namely that strength is conversely proportional to the raster angle.
The two numerical methodologies demonstrated in this study depict specific advantages. The RVE-CLT method requires less work in numerical software, as it contains only 6 simulations, whereas the angular RVE needed 36 × 6 = 216 simulations. On the other hand, the latter omits the issue with the shear modulus that arose in the RVE-CLT method. The Finite Element Method is based strictly on the accuracy of the discretization of a 3D model, in this case, the RVE. A higher number of elements was implemented in the second approach; nevertheless, the simplicity of the geometry of the RVE proved that even in the RVE-CLT method, with fewer elements, the influence of discretization on the final results can be omitted.
The highest yet acceptable discrepancy among values of Young's moduli, equal to 12.29%, was demonstrated for the RVE-CLT approach. Numerical results that were more fitting to experimental procedures were achieved with the angular RVE approach, with the highest relative error just below 8%. Conducting several simulations with a parametrized angle of force applied increased the accuracy of the results. The smallest errors were noted in the S1 and S2 series, i.e., without interchanging layers, where experimental S2, compared to both numerical methods, depicted a higher error equal to 6.38%.
Data acquired from the experimental test was used as input to an algorithm based on Classical Lamination Theory. The direct approach did not provide accurate results only in the case of the shear modulus, maintaining accuracy for Young's moduli and Poisson's ratios. Great attention has to be paid while analyzing analytical results, as could have been seen in this case-if not for Saint-Venant's Principle and redefining the equation for shear modulus, the results would be inaccurate, especially since the initial stress was smaller.
Regarding the incorrect breaking of the specimens, selecting a more suitable mechanical test that would still provide necessary mechanical properties, such as inter-laminar shear or three-point bending, could overcome the problem. Analytical as well as numerical methodology ought to be adjusted in this case.
One of the most crucial factors influencing overall strength in any case investigated in this study is the extent of the void inside the sample. Fibres in the laminate, in this case, filament strips, did not fully overlap, which caused the formation of gaps between them. Decreasing the gaps' size could positively influence the material's mechanical properties, i.e., tensile strength or Young's modulus. Melting temperature or general printing speed during the specimens' manufacturing could have influenced the microstructure of the samples. For instance, a lower speed could improve gap formation. Nevertheless, fully removing gaps in FFF 3D printing brings various issues to the printing process. The microstructure shows that the cross-section of a filament strip is not a regular polygon; hence, it will not fully tessellate without creating gaps. In order to build a microstructure without gaps, the strips ought to overlap significantly, which, in turn, vastly increases printing time.
When it comes to measurement reading, the tensile test machine used in this study provided data that had to be analyzed and computed to obtain mechanical properties. Even though the machine testing range was 50 kN, offering lower accuracy in the force interval under inspection, the differences in behavior between samples with different raster angles were well visible. Later, proper mathematical algorithms prepared in the programming language Python diminished the labor work connected with the determination of these properties. On the other hand, technologies such as DIC (Digital Image Correlation) can provide an entire map of the strain field within the tested structure without direct contact with the sample, compared to the extensometer used in this study. Such reading could offer better conclusions on how strains are displaced within the sample while undergoing tensile stress.
Fused Filament Fabrication is becoming increasingly popular, even among private individuals. Since most do not have an advanced understanding of micromechanics, a rule was firmly established: G-codes are keen to be created with a 45 • filament orientation. If such a model undergoes any stress, it is not pure tension, as it was tested in this study; nevertheless, it still showed that [45,−45] have a specific characteristic, that is, the highest strain before breaking.

Conclusions
In this study, 3D printed structures are tested experimentally in terms of their response to tensile stress, and numerical and analytical approaches are made to predict strength behavior in a spectrum of raster angles of filament strips in the structures. Experimental measurements-the tensile test chosen for this study was conducted, attaining accurate results but with incorrect breaking. The numerical part of the study proved that direct results taken from FEM software were unreliable and had to be altered through a nonsensitive shear modulus methodology. In addition, a more profound simulation was carried out, including parametrization. Finally, the results from all procedures were compared and analyzed. The RVE-CLT approach showed a maximum discrepancy in Young's modulus values of 12.29%. A closer match to experimental outcomes was found using the angular RVE method, with errors up to 8% as the method turned out to be non-sensitive to the determination of shear's moduli. Also, the study has shown that material behavior greatly depends on building parameters, mainly filament orientation. Contrary to the isotropic material type of PLA, 3D-printed specimens made with it exhibited strictly orthotropic behavior.
Research conducted in this study contributes to the determination of a suitable approach to simulate stress behavior in FFF 3D printed structures numerically. Such studies are highly important, considering the development of additive manufacturing technologies and other technologies, such as SLA or SLS. Even though the main goal of this study was achieved, there is a wide range of future opportunities for this study, such as using more sophisticated material characterization techniques (e.g., scanning electron microscope, hole-drilling technique, or dynamic mechanical analysis, allowing for mesoscale modeling). Failure simulation approaches could simulate breakage. As mentioned earlier, the [45/−45] layup exhibited particularly high strain at break that could be studied more profoundly with other layups, especially those near 45 • or −45 • .

Data Availability Statement:
The data that support the findings of this study are available from the corresponding authors upon request.