1 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.

Fig. 1
figure 1

Drawing of simplified cellular material structure of 4 filament lines placed alongside another, creating a void in the enclosed center region. Longitudinal and transverse directions, l and t, are defined relative to the axis of printing

Full size image

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).

2 Voids in fused deposition modeling

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 Fig. 2:

  • Reduction of cross section due to voids;

  • Void-induced stress concentrations;

  • Incomplete interdiffusion of polymer chains.

Fig. 2
figure 2

Strength reduction mechanisms for transverse loading in FDM manufactured specimens

Full size image

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.

Fig. 3
figure 3

(a) Virgin sample, (b) failure due to longitudinal loading (filament fracture), (c) failure due to transverse loading (inter-filament fracture), and (d) failure due to vertical loading (inter-layer fracture).

Full size image

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 wash-boarding effect. This is especially evident for first-layer 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.

Fig. 4
figure 4

First-layer-defects due to (a) excessive accumulation of material where the nozzle is too close to the print bed, generating a wash-boarding effect, and (b) insufficient infill due to too large distance between nozzle and bed

Full size image

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).

Fig. 5
figure 5

Shape of voids from CT scan data, with an assumed model of the filament boundaries

Full size image

3 Method of CT data analysis

The following procedure extracts information about the sizes of voids in a sample of cellular material where the voids are:

  1. 1.

    Known in number.

  2. 2.

    Periodically distributed (as opposed to randomly distributed).

  3. 3.

    Possible to isolate (each void is contained by four adjacent lines of material).

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,

Fig. 6
figure 6

Procedure for measuring voids: (1) Obtain volumetric data from specimen. (2) Threshold the data/images, separating air and solid material. (3) Set up a cell-grid, where each cell covers an area where one would expect one void. (4) Characterize the height and width of each void. (5) Identify the largest void in each cell. (6) Analyze the through-thickness properties

Full size image

applies the following procedure:

  1. 1.

    Obtain volumetric data from specimen.

  2. 2.

    Threshold the data/images, separating air and solid material.

  3. 3.

    Set up a cell-grid, where each cell covers an individual intersection between four lines of filament

  4. 4.

    Characterize the height and width of each void.

  5. 5.

    Identify the largest void in each cell.

  6. 6.

    Analyze the through-thickness properties.

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:

$$ \left[\mathrm{x},\mathrm{z}\right]={\displaystyle \begin{array}{c}+{\mathrm{N}}_1\left(\upeta -1\right)\left(\upxi -1\right)\left(\uppsi -1\right)\\ {}+{\mathrm{N}}_2\left(\upeta \right)\left(\upxi -1\right)\left(\uppsi -1\right)\\ {}+{\mathrm{N}}_3\left(\upeta \right)\left(\upxi \right)\left(\uppsi -1\right)\\ {}+{\mathrm{N}}_4\left(\upeta -1\right)\left(\upxi \right)\left(\uppsi -1\right)\\ {}+{\mathrm{N}}_5\left(\upeta -1\right)\left(\upxi -1\right)\left(\uppsi \right)\\ {}+{\mathrm{N}}_6\left(\upeta \right)\left(\upxi -1\right)\left(\uppsi \right)\\ {}+{\mathrm{N}}_7\left(\upeta \right)\left(\upxi \right)\left(\uppsi \right)\\ {}+{\mathrm{N}}_8\left(\upeta -1\right)\left(\upxi \right)\left(\uppsi \right)\end{array}} $$
(1)

where the void row, columns, and image number relate to these values in the following ways:

$$ \eta =\frac{i-1}{l-1} $$
(2)
$$ \xi =\frac{k-1}{m-1} $$
(3)
$$ \psi =\frac{j-n}{p-n\ } $$
(4)
Fig. 7
figure 7

Linear interpolation of coordinates for setting up a grid for finding voids. The nodes N are the [x, z] coordinates of the outermost voids in the two slices used for interpolation

Full size image

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.

Fig. 8
figure 8

The column, slice, and row index of the voids

Full size image

As seen in Fig. 9, for all calculated void positions, a slice of the dataset with height and width, corresponding to 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.

Fig. 9
figure 9

CT slice with indications for each interpolated void position, and a rectangular box indicating the current cell, together with detail of a cell and its void maxima

Full size image

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:

$$ {C}_h^{ijk}=\frac{h}{H}=\frac{z_4-{z}_2}{H} $$
(5)
$$ {C}_w^{ijk}=\frac{w}{W}=\frac{x_3-{x}_1}{W} $$
(6)

where [xn, zn] are the coordinates from points Pn for the voids in cell ijk. The area was simplified, using the outer points:

$$ {a}^{ijk}=\frac{x_1{z}_1-{x}_2{z}_1+{x}_2{z}_3-{x}_3{z}_2+{x}_3{z}_4-{x}_4z+{x}_4{z}_1-{x}_1z}{2} $$
(7)
$$ {C}_a^{ijk}=\frac{a^{ijk}}{HW} $$
(8)

When written as a one or two index average value, e.g., \( {\overline{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.

Fig. 10
figure 10

Position of start and stop frame/slice for analysis, as seen from a horizontal slice of the specimen

Full size image

4 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 approximate relationship between the rotational velocity of the feeder and volumetric extrusion rate of material, experimental results from Bellini et al. [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.

5 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. 1.

    Use the previously described method to obtain the shape and area of the inherent voids.

  2. 2.

    Find the distribution of void sizes and average void shape.

  3. 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. 4.

    Map these onto a finite element model using the distribution data from step 1.

  5. 5.

    Simulate the uniaxial loading response of the model.

  6. 6.

    Find the global stress/strain response and energy storage in each element.

Fig. 11
figure 11

Proposed method for finding stiffness of a uniaxially printed FDM part, and its internal stresses consisting of (1) processing of CT-data; (2) generating characteristic void shapes; (3) finding unit cell responses; (4) mapping unit cell data onto global structure; (5) simulating stresses and strains (6) calculating global response and strain energy density distribution

Full size image

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 the q-th corner of the void with index \( ijk,\kern0.5em \left[{x}_q^{ijk},{z}_q^{ijk}\right] \), the average relative position \( \left[{\overline{\Delta x}}_q,\overline{\Delta {z}_q}\right] \) using \( \left[{x}_1^{ijk},{z}_1^{ijk}\right] \) as reference point in each void (excluding cells without voids) is calculated as:

$$ \left[{\overline{\varDelta x}}_q,{\overline{\varDelta}}_q\right]=\underset{i,j,k=1,1,1}{\overset{l,m,n}{\Sigma}}\frac{\left[{x}_q^{ijk},{z}_q^{ijk}\right]-\left[{x}_1^{ijk},{z}_1^{ijk}\right]}{ijk-{N}_{\mathrm{zero}}} $$
(9)

where Nzero 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 Ca 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:

$$ \frac{\overline{\varDelta {x}_1}+\overline{\varDelta {x}_3}}{2}=\frac{W}{2} $$
(10)
$$ \frac{\overline{\varDelta {z}_2}+\overline{\varDelta {z}_4}}{2}=\frac{H}{2} $$
(11)

The unit cell response was found through a first-order 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:

$$ {u}_i^A-{u}_i^B={\varepsilon}_{ji}\left({x}_j^A-{x}_j^B\right) $$
(12)

where xi is the i-coordinate of the surface, ui is the displacement of the surface in the i-direction and εji are the imposed macro-strains applied as six different deformation modes.

$$ {\varepsilon}_{1\dots 6}=\left\{\left[1\ 0\ 0\ 0\ 0\ 0\right],\dots, \left[0\ 0\ 0\ 0\ 0\ 1\right]\right\} $$
(13)

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:

$$ {\overline{\sigma}}_{ij}=\frac{1}{V}\int {\sigma}_{ij} dV=\frac{1}{V}\sum \limits_{n=1}^N{\sigma}_{ij,n}{V}_n $$
(14)

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 Ca 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.

Fig. 12
figure 12

Loads/boundary conditions for loading in x-direction

Full size image

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 approaches using strain energy release rates as basis for fracture in linear elastic materials, while Rice [20] developed approaches for non-linear elastic and elastic-plastic 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, unorm for each cell i, j, k, could be found through:

$$ {u}_{norm}^{ijk}=\frac{u^{ijk}}{\frac{1}{2}\ F{\delta}_{bulk}} $$
(15)

where uijk 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:

$$ {\delta}_{\mathrm{bulk}}=F\frac{l_n}{E{A}_n} $$
(16)

with Young’s modulus E, length and cross section of cube in loading direction (n), ln and An, respectively. Stiffness reductions factors R are calculated by:

$$ R=1-\frac{F\ {l}_n}{E{A}_n{\delta}_u} $$
(17)

where δu is the unit displacement applied in the simulation.

6 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/s2. 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.

Fig. 13
figure 13

Drawing of the printed cubes used in this study. A small pedestal is added to give some spacing between the specimens and the rotating base while CT-scanning. Specimens are scanned up-side-down rotating around the vertical axis. All values in mm

Full size image
Fig. 14
figure 14

Toolpaths for printing a single layer of the cubes. All layers start in the lower left corner. All specimens were placed at the center of the print bed

Full size image
Table 1 Print variations for the three samples
Full size table

7 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.

Fig. 15
figure 15

Detail from samples 1 and 2 showing large difference in void size

Full size image

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.

Fig. 16
figure 16

Through-specimen average void heights (\( {\overline{C}}_h^{ik} \), \( {\overline{C}}_h^{jk} \) and \( {\overline{C}}_h^{ij} \) for left, right and upper faces respectively) in percentage of layer heights. Specimens 1 to 3, from left to right. Raster pattern shown to the left

Full size image
Fig. 17
figure 17

Through-specimen average of void width as percentage of filament widths (\( {\overline{C}}_w^{\kern0.5em ik} \), \( {\overline{C}}_w^{\kern0.5em jk} \) and \( {\overline{C}}_w^{\kern0.5em ij} \) for left, right, and upper plane respectively)

Full size image
Fig. 18
figure 18

Through-specimen average of void area divided by unit cell area (\( {\overline{C}}_a^{\kern0.5em ik} \), \( {\overline{C}}_a^{\kern0.5em jk} \) and \( {\overline{C}}_a^{\kern0.5em ij} \) for left, right and upper plane respectively)

Full size image
Fig. 19
figure 19

x-y plane with \( {\overline{C}}_a^{ij} \) of specimen no. 1, start in upper left corner, and color scale as in 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

Full size image
Fig. 20
figure 20

\( {\overline{C}}_a^{ij} \) for voids with left and right turning point, displaying largest void sizes mid to end of void, opposite to turning point (previous figure split into odd and even column numbers)

Full size image
Fig. 21
figure 21

\( {\overline{C}}_a^j \) for voids with left (blue) and right (green) turning point (averaged values from previous illustration)

Full size image
Fig. 22
figure 22

Defects in end-of-print location for each layer

Full size image
Fig. 23
figure 23

\( {\overline{C}}_a \) illustrated using different color scale for sample 1 and 3 for higher resolution. Detail of fluctuating void sizes, possibly due to machine dynamics

Full size image
Fig. 24
figure 24

Average cross section reduction in the transverse (\( {\overline{C}}_h^{\kern0.5em i} \)), longitudinal (\( {\overline{C}}_a^{\kern0.5em j} \)), and vertical (\( {\overline{C}}_w^{\kern0.5em k} \)) direction for the three specimens (1 to 3). All positional values taken starts at origin in the coordinate system in Fig. 7, and are relative to specimen size

Full size image

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 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 \( {\overline{C}}_w^{\kern0.5em k} \)). 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 are indications of machine or extrusion dynamics playing a role, as there is some oscillating variation along each void, where \( {\overline{\mathrm{C}}}_{\mathrm{a}}^{ij} \) for the two densest samples show oscillating values of 2–3%, as seen in 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.

Table 2 Maximum reduction in cross section along each axis. Maximum values for transverse (\( {\overline{C}}_h^{\kern0.5em i} \)), longitudinal (\( {\overline{C}}_a^{\kern0.5em k} \)), and vertical (\( {\overline{C}}_w^{\kern0.5em j} \)) void measures
Full size table
Table 3 Total average coefficients
Full size table

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 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 Ca would be equal to the porosity.

8 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.

Fig. 25
figure 25

Average void shape, based on relative distances between the corner voids, excluding cells without voids. Dimensions are in pixels, at approximately 119 pixels per mm

Full size image
Fig. 26
figure 26

Histogram of the Ca values for sample 1

Full size image

Unit cell simulation models were made for 100 different void sizes, with Ca ranging from 0 to 0.16 (equivalent to void cross section areas from 0 to 0.0192 mm2) 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 z- and x-direction, as these are dependent on the void height and width Table 4.

Table 4 Stiffness reduction in each direction compared with bulk material
Full size table

The reduction in stiffness in z-direction, 3.8%, is the same as the value for \( {\overline{C}}_a \) for this specimen. For a geometry with uniform voids this should be expected as Ca 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.

Fig. 27
figure 27

Normalized strain energy density for loading in transverse direction (exterior only)

Full size image
Fig. 28
figure 28

Normalized strain energy density for loading in transverse direction. Each cross section (column index i) shown. Images are divided into odd and even numbers for illustrative purposes due to alternating compliant and stiff side

Full size image
Fig. 29
figure 29

Normalized strain energy density for loading in longitudinal direction (exterior only)

Full size image
Fig. 30
figure 30

Normalized strain energy density for loading in longitudinal direction. Cross section (slice index j) shown

Full size image
Fig. 31
figure 31

Normalized strain energy density for loading in vertical direction (exterior only)

Full size image
Fig. 32
figure 32

Normalized strain energy density for loading in vertical direction. Each cross section (row index k) shown

Full size image

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, x- and 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.

9 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.