A new method for assessing anisotropy in fused deposition modeled parts using computed tomography data

Voids in fused deposition modeled (FDM) parts are assumed to be a key driver for their anisotropic behavior. However, these assumptions are based on investigations of voids using only 2D data (microscopy images). This paper presents a new method to measure such voids by analyzing 3D-data of from X-ray computed tomography (CT), and application of this data for assessment of mechanical parameters. The article is divided into three parts, where the first part elaborates on a proposed method to assess and characterize the void geometry throughout uniaxial printed FDM parts using CT-data. The second part presents an investigation of the void configurations in samples manufactured using different process parameters, aiming to understand how variation in extrusion rate and compensation for non-linear dynamic extrusion behavior affects the void sizes. The third part displays how the information regarding void sizes could be further related to global mechanical properties, using a multiscale finite element approach. The present method of CT-data analysis gives a clear overview of the spatial variation of the void geometry, and findings from the investigated samples suggest that the size of voids have a large non-random spatial variation, highly dependent on the turning points of the toolpaths, and also significantly affected by accumulation of excess material. Printing at a low extrusion rate increases the void sizes considerably, while implementation of an extrusion dynamics compensation algorithm was found to have low impact on the void sizes. The multiscale finite element approach predicts anisotropic elastic behavior, significantly more compliant in the vertical and transversal direction, relative to the printing direction of the infill. It also predicts a non-isotropic strain energy density throughout the specimen, where the location and magnitude of the most energy dense locations vary significantly for different directions of loading, which implicates an anisotropic behavior in terms of failure, in accordance with literature.


Introduction
Components manufactured by FDM are found to result in an anisotropic cellular material behavior [1][2][3]. For low relative density components, the internal structure, also called infill, typically consists of either rectilinear webs, as variations of a (0°, 90°) or (0°, 60°, 120°) alternating raster structures, or honeycombs. For structures with relative density above 75-80%, the most common methods are rectilinear or concentric infills, and as the relative density approaches unity, the voids left in the structure appear mainly in terms of small channels made by aligning four filaments alongside each other as seen in Fig. 1.
The parameters governing the geometry of these voids are numerous, including anything from machine resolution and machine dynamics, to temperature-dependent transient material flow and toolpath placement. Attempts have been made in modeling some of the processes involved in fusing filament segments [4][5][6][7], but there is still a long way to go to model the whole system. However, it is possible to investigate the end result of the printed component. This has previously been done by both a scanning electron microscope (SEM) and microscopy [1,8]. Unfortunately, these methods do only capture the void sizes along a single surface, and through thickness variations are therefore usually neglected in further analysis. In a previous study, the statistical distribution of void sizes combined with a weakest link approach was found to be the main driver for anisotropic strength characteristics of tensile specimens in PLA [3]. The method used a microscopy picture of a center slice of a specimen to estimate the void sizes, and the distribution was then analyzed and the expected weakest link of a given sample length was predicted. The theoretical results corresponded well with experimental ones, yet the estimated failure load distribution was lower than that found in experimental results. One of the possible explanations for this under-prediction was that the voids in the center slice are not representative for the full spatial distribution of voids, hence leading to a non-accurate prediction of strength. Investigating this hypothesis is thus one of the main motives for conducting this research. The contributions of this paper can be stated as follows: & A method to assess void geometries throughout the specimen; & A demonstration of how this method can be used to assess influence of process parameters on the voids, and investigate the effect of altering flow rate and pressure advance parameters; & A demonstration of how this method can be used to assess the voids influence on mechanical parameters as stiffness and failure parameters.
The method presented analyzes volumetric data from X-ray computed tomography (CT-scan), and measures the void sizes throughout a small cubic volume of near-dense FDM parts. The data is analyzed using the OpenCV image processing library through Python programming. Volumetric data from parts made by FDM has been researched earlier with a focus on dimensional accuracy and porosity [9,10]. However, to our knowledge, this has not been done with the focus on voids as a driver for anisotropy. In response to this gap in prior research, this article reports on the first attempt in relating 3D void geometries based on experimental volumetric information, to anisotropy properties.
To demonstrate how this analysis can be used for obtaining mechanical parameters of the resulting cellular material, a multiscale simulation method is developed. This method is capable of capturing the global elastic behavior and local strain energy density distribution, which can further be used to explore mechanical aspects of failure in FDM manufactured specimens. This is achieved through a first-order homogenization method using finite element analysis, using the geometrical configuration of the voids found through the CT-data analysis.
The reminder of the article is structured as follows: Section 2 presents previous literature on void formation in FDM. Section 3 describes the invented method used for analysis. Section 4 gives a brief overview of the process parameters investigated. The multiscale modeling approach is presented in Section 5. Section 6 displays the attributes and procedure for manufacturing of the samples. Experimental results and discussion from void analysis are given in Section 7, while experimental results and discussion from multiscale analysis are presented in Section 8. Finally, summary and key takeaways are found in Section 9.
It is important to note that there is a large variation in the design and performance of FDM systems, so that numerical results would most possibly be different from machine to machines. But we expect the trends to be similar, especially for the most popular Open Source systems, building on the Marlin Firmware, and even more so if also using a direct extruder design with both hot-end and extruder co-located on the print head (as opposed to a Bowden design).  The most conventional method of manufacturing internal structures by FDM, stacking small strings of filament, creates a cellular material with unit cells in the meso-scale. The resulting cellular material properties are usually anisotropic [1,11]. There are three important strength reduction mechanisms for parts produced by FDM, illustrated in The reduction in solid cross section explains most of the discrepancy between transverse (extrusion paths are transverse to the loading direction) and longitudinal (extrusion paths along the loading direction) tensile capacity [3,12]. General findings show that transverse loading results in a plane inter-filament fracture, longitudinal loading creates a fracture surface transverse to the printing direction, while vertical loading separates the layers, as illustrated in Fig. 3 [1,8]. The bonding between layers is therefore the key governing mechanism determining tensile capacity.
In contrast to voids in cast/molded components, the occurrence of voids in parts made by FDM is of a highly repetitive and structured nature, as they are dictated by the toolpaths of the printer. However, the size and shape of these voids are influenced by several process parameters, as well as inaccuracies and variations inherited by the FDM process. Moreover, machine dynamics, resolution and accuracy of the setup, and complex material flow involving molten plastics in a near-solid-phase would all affect the size and shape of the voids. Furthermore, when the filaments are placed side-to-side, small fluctuations in material flow or unevenness of the print surface will propagate from one filament line to the next, generating a washboarding effect. This is especially evident for firstlayer printing, where the print quality is very dependent on the evenness of the print surface. In this connection, an uneven surface would often lead to excessive deposition of material in some regions, while leaving other regions insufficiently infilled, as seen in Fig. 4.
For layers located further away from the print bed, the filament shape is more consistent as the deposition surface (the previously printed layers) has been leveled out. The top of each layer is smoother than the bottom, resulting in voids that resemble up-side-down kites or triangles as reported in [3] (Fig. 5).  The method is based on analyzing computed tomography (CT) data, but would work for any dataset representing a voxel-based 3D structure with material data that could be used to separate voids from solid material. The proposed method, as shown in Fig. 6, applies the following procedure: 1. Obtain volumetric data from specimen.  This approach, finding the void geometries rather than the conventional method of analyzing the neck geometry [6,12,13] (distance between the voids), allows for more computationally efficient analysis, as one do not need to find the nearest neighboring voids.
To be able to extract continuous data about every single void, a coordinate mapping system is used to estimate the position of the voids. This is done to compensate for specimen misalignment during the CT scan. Two CT scan slices at 1/3 and 2/3 width positions of the cubes, taken normal to the longitudinal direction of voids/filament lines, is used to set up this coordinate system. The coordinates of the corner voids of the two slices are named N1 to N4, and N5 to N8, respectively, as displayed in Fig. 7. The coordinates of each void are then found through linear interpolation, as follows: where the void row, columns, and image number relate to these values in the following ways: Here, i, j, and k are the void's column, slice, and row number, respectively, as seen in Fig. 8, l and m are the total column and row numbers, respectively, while n and p are the numbers of the first and second frame used for setting up the coordinate system. The reason for not using the first and last frames for setting up the coordinate system is due to geometry inconsistencies in those areas, experienced through trial and error.
As seen in Fig. 9, for all calculated void positions, a slice of the dataset with height and width, corresponding to  Fig. 8 The column, slice, and row index of the voids the layer height and filament width, was scanned for voids. This was achieved by using Python and the OpenCV library functions mean threshold and find contour. If more than one contour was found in one cell, only the largest one was stored. Noise was filtered out by using a minimum required contour length of seven pixels.
The maximum extension in x and y direction was found as seen in Fig. 9, and displayed as fractions of layer height H and filament width W of the void: where [x n , z n ] are the coordinates from points P n for the voids in cell ijk. The area was simplified, using the outer points: When written as a one or two index average value, e.g., C ij denotes the average over the missing indexes, in this case, the row index k. As some cells have more than one void (in the case of much noise or irregular void shapes), all voids in each cell are analyzed, and only the maximum values for these coefficients are used.
The geometrical discontinuities at the edges were neglected by defining the boundaries 1 mm from the edges, as illustrated in Fig. 10.

Parameter selected for investigation
Two parameters are selected for a brief investigation of their influence on the void size and spatial distribution, namely flow rate and pressure advance.
Flow rate, also called extrusion multiplier, refers to relative increased or decreased feeding velocity of raw material/filament. Although the rotational velocity of the filament feeder, and its radius will give an  [14] show that this deviation could be large. Tuning of this parameter is therefore crucial to achieve the correct amount of the extruded material [15]. This parameter  is highly connected with tensile capacity [1,11], and is therefore assumed to have a significant influence on the void sizes.
Pressure advance is a compensation algorithm for reducing unwanted extrusion defects found in areas of high acceleration or deacceleration and is therefore assumed to decrease eventual variation in void size near the edges. The specific algorithm used in this research was linear advance v1.0 as implemented in the Marlin v1.1.8 firmware. The algorithms' influence on the mesostructure or mechanical performance of FDM parts has not been investigated earlier.

Multiscale simulation approach
As the material is clearly non-uniform, the part's response to uniform (axial) loading is also expected be non-uniform. As a proposal of a multi-level approach for finding the global behavior, and generating local responses of strain energy density, we have adopted a multiscale finite element-based simulation strategy, as illustrated in Fig. 11. The approach includes the following steps: 1. Use the previously described method to obtain the shape and area of the inherent voids. 2. Find the distribution of void sizes and average void shape.
3. Create unit cell simulation models for cells with a range of void sizes, but same shape as found in step 2, and obtain effective homogenized material constants. 4. Map these onto a finite element model using the distribution data from step 1. 5. Simulate the uniaxial loading response of the model. 6. Find the global stress/strain response and energy storage in each element.
In theory, one could have made a simulation model for each void shape and size, but this would be computationally expensive. Therefore, to increase computational efficiency, a characteristic shape was used to represent the void geometry, found through averaging the relative positions of the corners for all voids. Naming the x and y position of   where N zero is the number of cells where the processing method did not identify any voids. Unit cell models were then developed by scaling this characteristic void shape according to the values for C a found through the CT data processing. It was assumed that the location of the voids inside each unit cell could be approximated by positioning the average of the x and y maximums in the center of the cell as follows: The unit cell response was found through a firstorder homogenization approach, as described by Geers et al. [16] using periodic boundary conditions and linear stress responses. The method assumes an insignificant strain gradient across the homogenized structure. Periodic boundary conditions for the unit cell response were implemented by identifying pairs of parallel surfaces {A, B} on the unit cell, and applying the following constraint equations: where x i is the i-coordinate of the surface, u i is the displacement of the surface in the i-direction and ε ji are the imposed macro-strains applied as six different deformation modes.
i and j should in this context not be confused with the column and image index previously described. The resulting average stresses from solutions of the individual load cases applied on the unit cell were calculated by multiplying the stresses at the centroid of each element by its volume fraction relative to the unit cell volume, and summing all elements in the model as shown: The finite element analysis was performed using the commercial finite element code ABAQUS, release 6.14. Models were created using the Python scripting interface supported in ABAQUS/CAE meshed with 8-node elements of type C3D8R. The resulting stress vectors represent the rows of the stiffness matrix, and the resulting effective elastic properties were extracted from the compliance matrix by inverting the stiffness matrix. Poisson's ratio in the simulation was set to 0.35 [17]. While irrelevant due to normalization of results, the Youngs modulus was set to 3000 MPa. The effect of element size was investigated in a few convergence studies for selected void sizes, and was found to have insignificant effect on the displacement results on the boundaries. Mapping of the cellular material properties from the cell-level simulations onto the global finite element model was done with reduced resolution in y-direction, so that the model could be simulated with 26 × 100 × 39 elements rather than 26 × 1100 × 39, which was the resolution from the CT scan data analysis. The C a used for each element was therefore the average over neighboring ± 5 data points taken in the y-direction, which was done to reduce the computational load.
The loading of the global simulation model was applied as unit displacement between the faces normal to the direction of loading, and these faces were also constrained from out-of-plane displacements. An illustration of the constraints for x-directional loading is shown in Fig. 12. The results were then used for obtaining the global stiffness by obtaining the resultant forces.
For relating the result to strength parameters through fracture mechanics, strain energy density is found to be a key factor. Griffith and Irwin [18,19] proposed  [20] developed approaches for non-linear elastic and elasticplastic behavior. Although this article will not go the full length and provide predictive approaches for failure in FDM parts, we will provide the strain energy distribution of the loaded structures. The strain energy density results are normalized, so that the strain energy density in each element is divided by the average strain energy density for an isotropic homogenous cube of with bulk material properties, loaded with the same resultant force. For loading in the n-direction, this normalized strain energy density, u norm for each cell i, j, k, could be found through: where u ijk is the strain energy density for each element, while F is the forces acting on the boundaries of the system and δ bulk is the boundary displacement for a dense cube of the same bulk material, loaded with same resultant force. This can be found through: with Young's modulus E, length and cross section of cube in loading direction (n), l n and A n , respectively. Stiffness reductions factors R are calculated by: where δ u is the unit displacement applied in the simulation.

Manufacturing and CT scanning of samples
The cubes used in this analysis were produced with dimensions as seen in Fig. 13, which is a trade-off between capturing a great number of voids and achieving a decent spatial resolution from the CT reconstruction. The cubes are printed with the same uniaxial toolpaths for each layer, as shown in Fig. 14. Three samples were printed in transparent polylactic acid filament provided by Flashforge, with variation in flow rate F, and pressure advance parameter K as seen in Table 1. Printing temperature was set to 210°C, print speed of 48 mm/s, and an acceleration of 2000 mm/s 2 . All samples were printed on a Prusa i3 MK2S, and due to firmware conventions, the pressure advance parameter is applied with a scaling factor of 512. All samples were scanned using a NIKON XT H 225 ST X-ray CT scanner with 1571 projections.

Void analysis results and discussion
By visual reference only, there is clearly a large difference on the void sizes between samples with flow rate 1 and 0.9, as seen in Fig. 15, while there is little difference when including pressure advance. The analysis results comparing the different specimens are seen in Figs. 16, 17, 18, 19, 20, 21, 22, 23, and 24. For the three-dimensional illustrations in Figs. 16,17,18,19,20,21,22, and 23, the faces of the cubes are showing the through-thickness average values in this direction.
The first observation is that the void sizes are not randomly distributed, as they show clear signs of spatial dependency through non-random patterns both in transverse, longitudinal, and vertical directions. The  Fig. 20. Results show alternating side of toolpath direction-change affecting void sizes, which creates an alternating strong bond (+) and weak bond (−) on opposite sides of the specimen non-random size distribution also implies non-uniform stress distribution during uniaxial loading-especially for the force flow for transverse loading. The results also show that the voids are smaller close to the print bed than further up, especially for the 3-5 first layers, which is in accordance with literature [6] (also illustrated in Fig. 24 as C k w ). As seen from Fig. 19, the void size at the sides where the toolpath changes direction is significantly higher than for the rest of the structure, and creates a relatively stiffer bond on those sides. As the boundary of the analysis is 1 mm from the edge of the specimen, the actual solid turning point is not a part of the analysis. It is therefore believed that the temperature of the end of the previous line is sufficiently high, so that it sinters together with the new line to a higher degree than the more distant areas. In Fig. 19, the void size at the sides where the toolpath changes direction is significantly higher than for the rest of the structure, creating a relatively stiffer bond on those sides. As the boundary of the analysis is 1 mm from the edge of the specimen, the actual solid turning point is not a part of the analysis. It is therefore assumed that the temperature of the end of the previous line is sufficiently high, so that it sinters together with the new line to a higher degree than in the more distant areas.
As seen in Figs. 19, 20, and 21, the voids are increasing in size from the turning point until approximately the mid-plane of the specimen where it reaches a plateau. The oscillation in void size is also significantly higher near the middle of the sample. It is suspected that this is due to the fact that the velocity of the printer is higher in these areas. There is also a high void size at the end-of-print for all specimens, which is assumed to be a print defect due to effects while stopping the material extrusion and removing the nozzle.
For the two specimens with default extruder multiplier, results from the x-z plane show that the void sizes decrease considerably throughout each layer. The most plausible explanation for this effect is accumulation of excess material due to inaccuracies. There  Fig. 23.
Average values for void height, area, and width for each cross section normal to longitudinal, transverse and vertical direction respectively is shown in Fig. 24, and total maximum and average values are shown in Tables 2 and 3.
Comparing samples 1 and 3 reveals that incorporating a linear pressure advance of 0.06 s does not impact the void sizes to a large extent, but might increase consistency along the longitudinal axis. Figure 23 shows that the samples incorporating this are slightly less dense in the turning points of the toolpaths (slightly less blue on the right-facing planes). The difference is however marginal, so no firm conclusions are made. Comparing samples 1 and 2, however, reveals that under-extrusion of 10% impacts the void sizes considerably. It triples the void's cross section, while increasing the height and width approximately 70%, compared  Fig. 7, and are relative to specimen size with sample 1. Tuning the flow rate correctly is therefore a crucial task for achieving the smallest voids possible, which has been emphasized in literature [11,21,22]. Total averages are seen in Table 3, where C a would be equal to the porosity.

Results and discussion of multiscale approach
The average void shape is shown in Fig. 25. A histogram of the void size values is shown in Fig. 26.
Unit cell simulation models were made for 100 different void sizes, with C a ranging from 0 to 0.16 (equivalent to void cross section areas from 0 to 0.0192 mm 2 ) where the void shape was a scaled version of the geometry seen in Fig. 25, uniformly scaled in x and z-direction. These were then mapped onto the global structure.
The results from the simulations show a reduction in stiffness in all directions, with the largest reduction in zand x-direction, as these are dependent on the void height and width Table 4. The reduction in stiffness in z-direction, 3.8%, is the same as the value for C a for this specimen. For a geometry with uniform voids this should be expected as C a would be equal to the reduction in cross section for that direction, and hence the reduction in stiffness should be equal [19]. The non-uniform void size that these parts exhibit does, however, seem to not affect this value.
Normalized strain energy densities for transverse (x), longitudinal (y), and vertical (z) loading are shown in Figs. 27, 28, 29, 30, 31, and 32, displaying a considerable spatial dependency. Note that there are different legend scales for the three load cases. Based on the assumptions of planar fracture surfaces as shown in Figs. 3, 28, 30, and 32 display the strain energy densities for each cross section perpendicular to the direction of loading.
From the results, it is very clear that one introduces more energy into the system for loading in transverse or vertical direction compared with loading in longitudinal, using the same force. While the most energy dense elements for y-directional loading has approximately 10% more energy than for loading an isotropic cube of bulk material, xand z-directional loading have 37% and 43%, respectively. This could explain much of the cellular material anisotropy reported from FDM specimens.
The regions of high strain energy density also differ from4 load case to load case. For transverse direction loading, the highest strain energy density levels are found on the compliant sides as discussed in Fig. 19.  The strain energy density is also highest for the first columns due to the large void sizes, and on the very last, because of the large defects during print. Conversely, for loading in vertical direction, the high strain energy density levels are found on the stiff sides. For loading in longitudinal direction, the high strain energy density elements are somewhat more randomly distributed, but tend to be higher in the stiff region in the end of each layer (right hand side for each image in Fig. 30).
An important aspect is that many of the high strain energy density areas are located at or near the boundaries of the geometry, especially for loading in the transverse and vertical directions. This observation indicates that edge effects might play a considerable role for eventual crack initiation and growth.

Conclusions and further work
This article has presented a novel method for automatically capturing void sizes from CT scan data. In accordance with literature, the method is capable of identifying effects of decreasing flow rate on the void sizes, which increases the void sizes significantly, from approximately 27 and 28% of the cell height and width, to 47 and 49% with a 10% decrease in flow rate. The results also indicate that incorporating pressure control using an advance algorithm increases the void sizes marginally, but not enough to firmly conclude on any relationship. The key takeaway from the results, regarding void sizes, is that they exhibit a clear spatial dependency of void sizes; hence, prior research must be approached with care as most often single microscopy images are used for geometry assessments, and through-thickness variations are thereby omitted.
Further, the research has shown how geometry data of the voids can be used in investigating cellular material properties, linear elastic properties and strain energy density distribution through a novel multiscale approach based on first-order homogenization. This approach is applied on the default material sample under no pressure advance compensation or altering of flow rate. The results show that the introduction of voids make the structure significantly more compliant for loading in the vertical and transverse direction, resulting in 13.5% and 10.7% lower stiffness, respectively, than fully dense cubes with the same material. For longitudinal direction loading, the structure is 3.8% less stiff, which is the same magnitude as the porosity for the sample. The results from the multiscale model also show that the strain energy densities are much higher for loading in the vertical and transverse directions compared with loading in the longitudinal direction. The highest cell-average strain energy densities for the sample investigated is 37% (transverse loading), 11% (longitudinal loading), and 41% (vertical loading) higher than it would have been for a non-porous sample with the same loading. This indicates anisotropic failure properties, but a more complete failure assessment model has yet to be developed. Further work on the method of capturing the geometry of voids should try to optimize process parameters for  minimizing void sizes. Here, it is of interest to investigate whether this has any influence on the strength of FDM specimens by making comparison with experimental data. Another important extension of this work on void analysis would be to further develop the method to be able to analyze CT data from geometries with an alternating layup (e.g., [0°, 90°]) or more complex global geometries.  As failure in unidirectionally printed FDM, specimens would propagate in the interface between layers or columns of filament lines for transverse or vertical loading, further work on the multiscale simulation method should try to establish methods for identifying the weakest layer or column intersection. The approach could also be enhanced to include more geometric variation of the voids in the unit cell response simulation, rather than relying on a characteristic shape. However, as most of the high-strain energy density areas are on the boundaries of the geometry, a complete failure assessment should aim to include edge effects, as they could be essential for crack initiation When such a method is in place, each cross section can be analyzed individually, and be used to map the strain field from the global simulation onto the identified weakest cross section, before further investigating the link between FDM mesostructure and failure.

Funding information This research is supported by The Research
Council of Norway through BIA project 235410/O30, and is done in collaboration with AquaFence AS. We greatly acknowledge their support.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.