To study the effect of the overestimation in the thickness and the quality of approaches for correction, 16 diamond unit cells are generated. By shifting the diamond structure proposed by Huber et al. (2014) by a quarter of a unit cell length in all three coordinate directions (Soyarslan et al., 2017), four ligaments with complete nodes at both ends are positioned in the center of the RVE. These core ligaments are later analyzed with respect to their thickness distribution by different algorithms, as they remain unaffected by cuts at the boundary of the RVE.

In what follows, the investigation of the mechanical behavior is limited to macroscopic compression, which is commonly used in experiments (Jin et al., 2009; Huber et al., 2014; Hu et al., 2016; Liu and Jin, 2017). The resulting macroscopic properties are only valid for this loading direction. Due to the inherent anisotropy in the diamond structure, the mechanical response can be different for compression, tension, and shear. The elastic properties though can be considered isotropic in tension and compression, because elastic properties per definition reflect small deformations. Furthermore, because of the perfect symmetry of the unit cell in x, y, and z-direction, isotropy in these directions is naturally given as long as the loading is consistently either tension or compression. Thus, the stress-strain curve will show perfect agreement for small strains, whereas with increasing strain, the stress-strain curves for tensile loading tends to rise faster compared to the curves for compression loading. Under tensile loading, the ligaments tend to align in loading direction (see Sun et al., 2013) and are able to bear higher loads compared to compression loading, where the ligaments deform like an s-shape due to bending (Huber et al., 2014). Therefore, the yield strength is slightly larger in tension than in compression and the difference is more pronounced for thin ligaments, because they align more easily in tensile direction like fibers. These mechanisms are demonstrated for two example structures G₁₁ and G₁₄ in Supplementary Section 2.3. For the scope of this work it is sufficient to concentrate on compression, because errors in the ligament geometry will be reflected similarly in all mechanical properties and loading scenarios. In what follows, we will investigate the errors in the thickness determination depending on the algorithm that is used and their correction. To this end, we use diamond structure consisting of identical ligaments with well-known geometry. Because of this replication, the macroscopic behavior of the structure gives an indication about the response of a single ligament that is part of a more complex network.

Variable ligament shapes are incorporated in form of a continuous parabolic-spherical shape introduced by Richert and Huber (2018), see Figure 10 therein. To incorporate also asymmetric ligament shapes observed by Richert and Huber (2018), the ends are defined by two different radii r_{end, l} and r_{end, r} for the left and right junction, respectively. The resulting gradient along the ligament with length l is included in Equation (1) through the parameter b. The locations x_{Q, l} and x_{Q, r} at which the parabolic shape transitions into the spherical parts of the ligament, are determined iteratively such that a smooth ligament with a tangential transition is achieved (see Richert and Huber, 2018).

The axial coordinate x has its origin in the mid of the ligament, such that the ligament mid radius is given by r_mid = c. For the in-depth study of the thickness determination and correction as well as their effect on the mechanical properties, the ligament shape is kept symmetric by setting b = 0. In this case, r_end = r_{end, l} = r_{end, r} and x_{Q, l} = − x_{Q, r}.

In what follows, the unit cell size a_UC is set to 1, i.e., all absolute lengths are given as fraction of the unit cell size. The 16 geometries are chosen to cover ratios of ligament mid to end radius r_mid/r_end from 0.5 to 1.25 in increments of 0.25. This is the relevant range of ligament shapes as identified from a 3D tomography of a NPG sample (Richert and Huber, 2018). As the second geometry parameter, the end radius was varied from r_end = 0.1 to 0.175 in increments of 0.025. Through the combination of these two parameters a large range of solid fractions is covered that exceeds the typical range of NPG samples from very low (φ_min ≈ 0.1) to very large values (φ_max ≈ 0.5). Based on the two chosen parameters r_mid and r_end, the parameter c in Equation (1) can be determined following Richert and Huber (2018).

Reference FEM beam and solid models are generated for all geometries defined in Table 1. A detailed description of how the reference FEM beam is created is given in Supplementary Section 1. The solid unit cells are built using PCL scripting in MSC Patran 2017 and, after Boolean operation on all ligament and junction volumes, are meshed in a single meshing operation with C3D10 three dimensional 10-node quadratic tetrahedron elements for (Abaqus, 2014). The number of elements range from 9,445 to 38,279 for structures with lowest (G₁₁) and highest solid fraction (G₄₄), respectively, with average element sizes of 0.05. The solid fractions given in Table 1 are obtained from the FEM solid model in Abaqus via the history output VOL. Examples for the most filigree structures with r_end = 0.1 are shown in Figure 3. Due to the small ligament diameter, these structures will show the highest sensitivity with respect to effects of voxel resolution, discretization, and the accuracy of the algorithms applied to these data.

In addition to the solid models that serve as common reference for all mechanical properties, FEM beam models with 20 beam elements per ligament of type B31 [two-node shear flexible Timoshenko beams in space; (Abaqus, 2014)] are built using the code developed by Huber (2018). The code is modified for assigning a variable ligament shape to the beam elements in dependence of their position relative to the mid of the ligaments.

For the mechanical properties, a Young's modulus of E_s = 80 GPa, a Poisson's ratio of ν = 0.42, a yield strength of σ_{y, s} = 500 MPa, and a work-hardening rate of E_T = 1,000 MPa are chosen. These parameters represent the mechanical behavior of the ligaments in NPG reasonably well (Huber et al., 2014; Hu et al., 2016; Roschning and Huber, 2016; Huber, 2018).

The translation of the ligament shape given in Equation (1) for a single ligament into a physical meaningful radius distribution for the interconnected structure is described in detail in Supplementary Section 2. Through the intersection of three convex ligaments, the actual size of the nodal mass increases to the value R, which is defined by the triple point—the point where the surfaces of three ligaments intersect. This surface point is closest to the center of the nodal mass. Therefore, all reference FEM beam models are based on the radius for the biggest sphere R, that fits in the nodal area. The corresponding radii are computed as distance from the center of the junction to the surface in direction of the triple point, which is found at an angle of 70.53° relative to the ligament axis. The value R is assigned to all elements positioned between the ligament end, which is the center of the nodal mass, to the axial position of the triple point T. This approach avoids case sensitivity and allows to compare the results from different models. All geometric parameters for the structures defined in Table 1 are provided in Supplementary Section 4, Supplementary Table 1. Supplementary Figure 6A shows that there is only a moderate effect in the macroscopic Young's modulus. For most ligaments, the stiffening is below 10%. However, for the yield strength shown in Supplementary Figure 6B, the incorporation of R becomes relevant for cylindrical and convex shaped ligaments, for which a strength increase by up to 20% and 40%, respectively, is achieved.

For a finite model size, the choice of the boundary conditions can significantly influence the material response significantly. Miehe and Koch (2002) showed for shearing of a composite microstructure modeled with 2D solid elements that prescribed displacement boundary conditions lead to a stiffer response compared to periodic boundary conditions. Diebels and Steeb (2002) showed that boundary layers of rotations form under simple shear of a foam leading to a size effect. In this study, we investigate the effect of errors in ligament geometry on macroscopic properties and effects of boundary conditions should be avoided. Therefore, the chosen boundary conditions emulate an infinite periodic microstructure. Due to the perfect symmetry of the diamond structure, all simulations can be based on one unit cell with prescribed displacement and rotation boundary conditions, for details see Supplementary Section 1.2. For the FEM beam model, this approach is equivalent to periodic boundary conditions, while it significantly simplifies the meshing of a 3D FEM solid model.

The displacement boundary conditions impose the known deformation behavior of the structure on all surface nodes using ^*EQUATION in Abaqus. To this end, nodes on planes x = 0, y = 0, and z = 0 are set to zero displacement normal to the corresponding plane. Nodes in the planes at coordinate x = 1, y = 1, and z = 1 are set to remain in a plane that is controlled by a dummy node. All nodes on the mid planes are forced to move half the displacement of the corresponding nodes in the plane at coordinate 1. Finally, in the beam models, all rotational degrees of freedom are set to zero for all surface nodes. As no displacement boundary conditions are applied to the five internal junction nodes within the RVE, these nodes are allowed to move and rotate without any constraint. Nevertheless, they behave identically to the nodes at the boundaries, which have their rotational degrees of freedom fixed, and accomplish a full periodicity of the stress and deformation field results. This indicated the correctness of the chosen boundary conditions being equivalent to periodic boundary conditions. More details are given in Supplementary Section 1.2 (see Supplementary Figures 2, 3).

For elastic computations, a compression strain of 1% is applied on the dummy node of plane z = 1; for predicting elastic-plastic stress-strain behavior, the structure is compressed by 20% strain using large deformation theory (NLGEOM) with a start increment of 0.001. The Young's modulus is always determined from the first loading increment.

For geometries G₁₁ and G₁₄ (r_end = 0.1, r_mid/r_end = 0.5 and r_mid/r_end = 1.25, respectively), a size study with RVEs of increasing model size confirmed that the chosen displacement boundary conditions yield results identical to periodic boundary conditions, both being independent of the model size. The results are presented in Supplementary Section 1.2. As shown in Supplementary Figure 3, the computations with simple symmetry conditions, as used e.g., by Huber et al. (2014), asymptotically approach this value with increasing model size (see also the size study in the Appendix of Huber, 2018). For applying the displacement boundary conditions in the solid model, a search tolerance of 1% of the unit cell allows collecting enough FE nodes, which are sufficiently close to the position of the corresponding surface nodes of the FEM beam model.

Figure 4 shows contour plots for the corresponding FEM solid and beam models at a deformation in the elastic-plastic transition. Elements exceeding the yield stress of 500 MPa are colored in gray. They represent the distribution of the plastic zones, which are in good agreement for the solid model and the corresponding beam model for the convex ligament shape G₁₄, as can be seen from Figures 4C,D. However, for structure G₁₁ with concave ligaments shown in Figure 4A, the plastic zones are organized in the FEM solid model along the tension and compression side in the thin regions of the ligaments and cross the junction volume in the middle into the neighboring ligament. Due to the kinematics implemented in the FE beam elements, the FEM beam model in Figure 4B cannot capture this complex deformation and localizes the plastic strains in elements in the transition region from the ligament to the nodal mass.

In the following section, the results obtained from the FEM beam model and the FEM solid model are presented for the reference geometries defined in Table 1. This serves two goals. The first goal is to precisely determine the differences between the macroscopic properties of the FEM beam model relative to the FEM solid model of the very same geometry for all ligament geometries. For all further investigations, the FEM beam models serve as reference for the FEM skeleton beam models derived from the voxel models. This allows to clearly separate potential effects from different sources, such as the different behavior of FEM beam and solid models, the thickness algorithms (section FEM Skeleton Beam Models), and the quality assessment of the developed correction methods (section Methods for Thickness Correction).

The macroscopic properties Young's modulus E and the yield strength σ_y are derived from engineering stress and strain measures (see Supplementary Section 1.2, subsection Macroscopic Evaluation). Complete sets of the resulting mechanical properties for the structures defined in Table 1 are provided in form of absolute values in Supplementary Section 4, Supplementary Tables 2–4. An overview of the macroscopic mechanical properties predicted by the reference FEM beam model (E(ref), σy(ref)) normalized to the corresponding values of the reference FEM solid model (E, σ_y) is given in Figure 5. The shaded regions indicate solid fractions that are out of the range of NPG (Liu and Jin, 2017; Soyarslan et al., 2018a). It should be noted that a direct comparison with NPG samples via the solid fraction is not possible, because a significant percentage of solid fraction can exist in form of dangling ligaments, whereas our diamond structure is fully connected. Therefore, the larger range of solid fractions in this theoretical work can be useful for covering the relevant ligament shapes determined by Richert and Huber (2018).

Figure 5A confirms that the FEM beam model generally underpredicts the macroscopic Young's modulus relative to the solid model, which is due to the well-known effect from increased lever length (Huber et al., 2014; Roschning and Huber, 2016). The FEM beam model is more compliant compared to the solid model, because the full distance from the mid of the element to the ligament end, i.e., the half ligament length l/2, is available for bending deformation, independent of the ligament thickness. In contrast to this, the nodal mass in the solid model reduces the lever length available for bending of the ligament depending on the size of the nodal mass relative to the ligament radius. The node is stiffened-up and deformation is moved into the transition zone from the ligament to the nodal mass. For more details, we refer to Huber et al. (2014) and Roschning and Huber (2016). Jiao and Huber (2017b) carried out a study on the effect of the nodal mass for a ball-and-stick model and suggested a nodal corrected beam model to compensate for the softening in the beam model by adjusting the radii and the Young's modulus of the elements in the nodal region.

There is a clear trend toward the stiffness of the solid model for decreasing ratio r_mid/r_end, which goes along with a decreasing solid fraction. This means that the more concave the ligament is, the closer the macroscopic mechanical stiffness is to that of the FEM solid model. Therefore, concave ligaments require less nodal correction to raise the stiffness by about 30% (r_mid/r_end = 0.5) or 80% (r_mid/r_end = 0.75), while cylindrical and convex ligaments require an additional stiffening by more than a factor of 2. This disproves an application of a single “stiffness intensity factor” as proposed by Soyarslan et al. (2018b) independent of the local ligament shape and solid fraction φ.

In contrast to the elastic behavior, the effect in the macroscopic strength, computed at 1% plastic strain, depends strongly on the specific ligament shape (see Figure 5B). In average, the yield strength predicted by the FEM beam model is comparable to that of the FEM solid model. However, for specific ligament shapes the ratio of the yield strength ranges from 0.6 to 1.6. An example is shown in Figures 4A,B. From the contour plots for both types of models it can be deduced that for concave ligament shapes, the plastic zone in the FEM solid model, Figure 4A, is distributed over a larger volume extending from one ligament via the nodal mass into the neighbor ligament. In contrast to this, for the FEM beam model shown in Figure 4B, the plastic deformation localizes in elements located in the transition zone from the ligament to the nodal mass. Therefore, the levers and resulting bending moments causing plastic deformation are longer in the solid model, effectively reducing its mechanical strength. This can explain the unexpected high strength of the FEM beam model for specific geometries.

Based on the good agreement of the yield strength averaged over all geometries, one could argue that a structure that contains a large range of ligament shapes does not require a nodal correction for the mechanical strength. This surprising result has important consequences for the interpretation of stress-strain curves predicted from FEM beam models derived from skeletonized structural data, because the elastic and plastic properties need to be treated differently.

To mimic the procedure according to the analysis of tomography data, a Python script is used to scan the reference RVEs for given ligament geometries. This scan produces a black (pore) and white (gold) tiff-stack in the chosen voxel resolution. Details on the tomography of the FEM beam models via parallel processing are provided in Supplementary Section 3. The tiff files of the 16 model geometries are available for download as the Data Sheet 2.zip folder of the Supplementary Material. Details of the files are provided in Supplementary Section 5. The code is validated using the open visualization tool Ovito by Stukowski (2010) confirming that the solid fraction of the voxelized model is below 1% error. To avoid boundary issues during the skeletonization and thickness analysis, as discussed by Richert and Huber (2018), a larger RVE of size 3 × 3 × 3 unit cells is used, similar to Soyarslan et al. (2018b). However, for the voxelization, the scan-box edge length around the mid-point is limited to 1.5 times of the unit cell size a_UC, so that on all sides exactly one additional ligament (0.25 of one unit cell) is connected to the center unit cell. The skeletonization is carried out on the resulting RVE of size 1.5 in the open-source software Fiji (Schindelin et al., 2012) with the BoneJ Plugin (Doube et al., 2010) Skeletonize 3D based on the thinning algorithm by Lee et al. (1994). The diameter estimation is carried out with the BoneJ Plugin Thickness (Dougherty and Kunzelmann, 2007) based on the Biggest Sphere algorithm by Hildebrand and Rüegsegger (1997) and the 3D Mathematical Morphology (TANGO) Plugin operation 3D Distance Transform by Ollion et al. (2013). The skeleton forms the beam element axis and the thickness data is used to calculate the section radii of the beam elements. For the FEM skeleton beam model building, only the data within the volume of the center unit cell is used. For further details about the procedure (see the Appendix of Richert and Huber, 2018).

The geometry G₁₁ with the smallest diameter was chosen to determine the accuracy as function of the voxel resolution. This most filigree structure with r_end = 0.1 and r_mid/r_end = 0.5 is shown in Figure 6A. Due to the small ligament diameter, it has the highest sensitivity with respect to effects of voxel resolution and beam discretization. The structure was scanned with 60, 100, 200, and 300 voxels per unit cell edge length N_v/a_UC (see Figure 6), yielding volume fractions of 9.2, 9.4, 8.0, 7.9%, showing a dependence on the voxel resolution. With the unit cell edge length a_UC = 1, one voxel has an edge length of 1/60 (0.0167), 1/100 (0.01), 1/200 (0.005), and 1/300 (0.0033) for the different resolutions, respectively. The smallest radius of the structure is 0.05 in the middle of the ligament. With the lowest resolution of 60 voxels per unit cell edge length, this results in only three voxels making up the ligament radius. With the resolution of 100 voxels shown in Figure 6B, the proposed minimum quality of 10 voxels per ligament diameter proposed by McCue et al. (2018) is met. The unsatisfying quality of the 60 voxels resolution leads to steps in the beam diameters and an uneven replication of the ligament profile, as visible in Figure 6A. As a consequence, local narrow neckings are averaged out, which leads to a stiffening of the mechanical response. In contrast, the 200 and 300 voxels resolutions show a satisfying quality of the surface (see Figures 6C,D).

When analyzing the effect of the different voxel resolutions on the mechanical behavior of the FEM skeleton beam models, the skeletonization, and originating from that, the discretization of the beam elements are further sources of errors. The skeleton of the structure is the one-voxel-wide centerline. It is achieved by surface thinning, as implemented in Fiji. Richert and Huber (2018) discuss different discretization approaches, where the most accurate approach appears to be to construct the beam axis as the connection between the centers of neighboring voxels (1 V/E). However, due to the discrete cubic size of a voxel, this can lead to harsh direction changes of up to 90° between two beam elements (zigzag). Especially for curved ligaments, as found in actual NPG tomography data, this has a great effect. This zigzag skeleton path results in a more compliant mechanical behavior, as shown by Richert and Huber (2018). The other approach is to average over a certain number of voxels. An approach of on average five voxels was tested by Richert and Huber (2018). This solves the issue with the skeleton zigzag on the one hand, but results in a lower number of beam elements per ligament on the other hand and, due to this, the ligament shape may be badly represented. Neckings are averaged out and the macroscopic stiffness and strength is probably overestimated. This is a similar effect as if using a low voxel resolution.

A new approach is introduced in this paper, were the skeleton voxels are fit by a Bezier function. This results in a smooth line, with the start- and end-node being fixed in their position. The Bezier fit is not forced to go exactly through the individual skeleton points of the ligament, so no overshoots arise, as is would be the case for a spline fit. The Bezier approach has the additional advantage that the desired number of equidistant beam elements per ligament can be chosen in dependent of the length and skeleton voxel number of the current ligament. For assuring comparability with the reference FEM beam model, 20 two-node shear flexible Timoshenko beam elements in space (B31) are used (see section Reference FEM Solid and Beam Models). For the boundary conditions (see section Boundary Conditions).

The models for the different discretization approaches based on the Thickness data are shown in Figure 7. The diamond structures analyzed in this paper have initially a straight ligament axis (Figure 7A). By using the discretization of one voxel per beam element (1 V/E) on a 60 voxels scanned structure, kinks are clearly visible in Figure 7B as tilted elements. Also for the 200 voxels scan resolution, the 1 V/E discretization shows kinks (see Figure 7C). This phenomenon is not avoidable due to the discrete voxel size, shape and orientation of the ligaments in space, even for the ideal geometries used in this work. This problem is solved via the newly introduced Bezier fit, which shows nicely aligned beam elements (see Figure 7D). Besides the discretization issues, the diameter overestimation through the Thickness algorithm is clearly visible in all three FEM skeleton beam models (Figures 7B–D), when compared to the reference geometry presented in Figure 7A.

The FEM skeleton beam model was built from the four different voxel resolutions of the geometry G₁₁ based on the Thickness diameter estimation algorithms. Furthermore, the two different discretization approaches with either each voxel being represented by one beam element (1 V/E), or a Bezier fit (Bez 20 E/L) are applied to the skeleton and diameter data. The results for the Young's modulus are displayed in Figure 8. The values are fit to a simple hyperbolic function E(Th)=k/(Nv/aUC)+E∞(Th), where the parameter E∞(Th) approximates the Young's modulus for a model with infinite number of voxels N_v/a_UC = ∞, as 590 MPa and 580 MPa for Bezier and 1 V/E discretization, respectively. The percentage deviation from those values is inscribed. The focus is here solely set on the effect of the voxel resolution and the two different beam element discretizations. The deviations to the reference beam model stemming from the diameter estimations are addressed in sections Thickness Analysis and Effect on Mechanical Properties.

Overall, the 1 VE models show slightly lower Young's modulus values than the Bezier models, and also lower deviations to its asymptotic value of 580 MPa at N_v/a_UC = ∞. As the skeleton is straight in the reference geometry, the effect of the increased compliance caused by the kinks in the ligament axis with the 1 VE discretization is small. For the lowest voxel resolution (60 V) the stiffness is overpredicted by up to 43% while for higher resolution, the accuracy increases. The Young's modulus of the models with a resolution of 200 voxels shows around 10% remaining difference to the predicted value at N_v/a_UC = ∞. Further refinement slowly increases the accuracy, but rapidly increases the computational time. Thus, all further computations will use the voxel resolution of 200 voxels per unit cell edge length N_v/a_UC with the Bezier fit to ensure comparability to the reference structures created with 20 elements per ligament. The remaining uncertainty in the prediction of the mechanical properties is up to 12% due to the voxel resolution and beam element discretization. The resulting voxel edge length of 1/200 (0.005) defines the achievable accuracy limit for the geometrical characterization in the following sections.

This section discusses the geometry derived with the Thickness algorithm (Th) and the Euclidean distance transformation (EDT) from the voxel scan of the underlying reference geometries, given in Table 1. Figure 9A shows the mean-radii 〈r^(.)〉 obtained from averaging over all 20 elements of a ligament normalized by the mean-radius of the reference geometry 〈r^(ref)〉. It can be seen that the deviation of 〈r^(.)〉/〈r^(ref)〉 increases with increasing concavity, independent of the end radius r_end and algorithm used. For the Thickness algorithm, the largest value of 1.2 is comparable to the results of Richert and Huber (2018), where values up to 1.3 have been reported using the mathematically exact ligament geometry as reference. It could be argued that the deviation of 20% in the geometry is still acceptable. However, as showed by Richert and Huber (2018), this causes serious overpredictions in the mechanical stiffness of the ligament by a factor of two. As expected, the data from the EDT show an underprediction for increasing concavity, but the relative deviations are significantly smaller compared to Thickness algorithm.

The advantage of the object-oriented-programming is that it enables to locally analyze parameters of individual ligaments at specific positions. Figures 9B–D show selected results for the effect on the local thickness determined in the end, quarter, and middle position of the ligament, respectively. From this series, the strength and weaknesses of each algorithm can discussed. In the overall comparison, the EDT algorithm is of superior accuracy. At the mid and end position, where the tangent of the ligament shape is flat, the diameter is determined with high accuracy. Only in the transition from end to mid position, represented by the quarter positions in Figure 9C, the expected underprediction can be seen in the EDT data. In the worst case that represents the largest diameter change, i.e., structure G₄₁, the deviation is −30%.

For the Thickness algorithm, the local overestimation of the r_mid value increasingly depends also on the absolute radius of the ligament end, the more concave the ligament is. This is a result of the following mechanism: The Thickness algorithm propagates the sizes of the nodal region into the ligament region. Firstly, all skeleton points inside the nodal sphere are assigned with this value R_node ≥ r_end, forming a nodal plateau of constant radius. Secondly, from there the ligament shape assumes a smooth transition from R_node to r_mid. However, in the extreme case of a very thick ligament, the two nodal spheres can even overlap in the mid position of the ligament. This would lead to an extension of the plateau over the whole ligament length. Due to this, the determined radius in the mid-point rmid(Th) can take all possible values from rmid(ref) to R_node.

In the following, we will investigate the impact of the determined geometries on the macroscopic mechanical properties. The question will be addressed in how far the averaged data or the local effects in the geometrical characterization are relevant in terms of the mechanical behavior.

In section Thickness Analysis, the deviations for the average and local thicknesses are determined for the 16 reference RVEs. Because the diameter enters the moment of inertia by a power of four in the stiffness calculation, the overestimation of the Young's modulus and yield strength is expected to be even higher. To quantify this effect, 16 FEM skeleton beam models are built from the 200 voxel resolution scans (section RVE Size and Voxelization), with a Bezier curve fit to the skeleton axis (section Skeletonization and Beam Discretization). In Figures 10A,B, the macroscopic Young's modulus E and the yield strength σ_y, respectively, obtained from the FEM skeleton beam model are compared to the values from the corresponding reference FEM beam model.

The factor of overestimation of the Young's modulus for the Thickness algorithm, presented in Figure 10A, is similar for structures with same ratio r_mid/r_end, independent of the absolute r_end value. Strongly concave structures show the highest deviations by up to a factor of 2. Tending toward cylindrical and convex structures, the deviation decreases to a factor of 1.2. The trend in the yield strength data in Figure 10B is similar, showing highest overestimations at strongly concave ligaments. With decreasing concavity, the decay is however more emphasized. Furthermore, stronger variations for different r_end values are observed, especially for the concave ligaments. There, smaller r_end values show higher overestimation, ranging from 1.68 to 2.15. The higher sensitivity of the yield strength is caused by the circumstance that the onset of plastic deformation results from the combination of weakest cross-section and applied bending moment, which again depends on the lever acting on this cross-section. In contrast to this, the elastic deformations spread over the whole ligaments and into the junction volumes and are therefore less sensitive to the local geometry (Huber et al., 2014).

It should be noted that cylindrical and convex ligaments show overall the lowest overestimation, which is still about 20% for both macroscopic properties. This is astonishing, as one might imagine that a cylindrical ligament should be perfectly reproduced by the Thickness algorithm. However, this is only true for a cylindrical ligament of infinite length. For the interconnected structure, which contains junction volumes that are larger than the cylindrical ligaments, the overestimation in mechanical properties is due to the mechanism discussed in section Thickness Analysis.

In line with the findings from the geometric analysis presented in Figure 9, the predicted deviations in the macroscopic mechanical properties for the EDT data are much smaller compared those obtained for the Thickness algorithm. The results can be considered accurate for cylindrical and convex shapes while for concave shapes the stiffness and strength are reduced up to 20%. If this is acceptable, the EDT can be used without further correction. It should be noted that stronger concavities or asymmetries as well as non-circular cross-sections can further increase these deviations also for the EDT.

Coming from the Biggest Sphere Thickness approach by Hildebrand and Rüegsegger (1997), the idea is to compensate its systematic trend of overestimation by the opposing equivalent, which is the Smallest Sphere approach, discussed by Liu et al. (2014). Between these two extremes, a Smallest Ellipse (SE) approach can be considered, as schematically presented in Figure 11. As input data, the coordinates of the medial axis and the respective Thickness values are used. Each point x along the axis located in a smallest ellipse inscribed into the Thickness data r^(Th), is assigned with the value of the ellipse major axis as the radius r^(SE). The ellipse allows to incorporate some flexibility in the range of radius assignment near the point under investigation. To this end, the linear eccentricity e of the ellipse was determined independently of the geometries defined in Table 1. Approximately 800 ligament geometries were created, reproducing the range of ligament geometries detected by Richert and Huber (2018), including asymmetric ligament shapes. A linear eccentricity of e = 0.75 produced the lowest errors.

The obvious drawback of the proposed Smallest Ellipse approach is that the minimum diameter of a ligament is bound to the minimum Biggest Sphere Thickness value. This can be seen in Figure 11, where in the center of the ligament a gap between the minimum radius of the original reference geometry and the reconstructed radius remains. In the nodal areas, the Biggest Sphere Thickness value represents the upper limit, which is correctly reproduced. The algorithm is robust since it does not require any assumption on a model function, parameter bounds, and parameter start values and works for symmetric and asymmetric ligament shapes.

The correction of the geometries via the Smallest Ellipse approach lead to an overall improvement in the predicted macroscopic mechanical properties (see Figure 12). The previously observed overestimation from 1.2 to 2.0 based on the Thickness data (see Figure 10A) is now reduced to an almost constant value between 1.1 and 1.25, i.e., the concave ligaments are most improved. As discussed before, the yield strength shows some stronger sensitivity to the different ligament shape parameters, while the overall improvement is comparable to that of the Young's modulus. In summary, the reconstruction of the ligament shape with the simple Smallest Ellipse approach represents a substantial improvement in comparison to the Thickness data, although some geometrical inaccuracies remain in the thinner region.

In contrast to the Smallest Ellipse approach, which does not require an assumption with regard to the ligament geometry, computational methods, such as optimization strategies or artificial neural networks can be applied for reconstruction of the original ligament shape. In our case, the ligament shape is limited to symmetric-parabolic shapes, which in principle allows applying both strategies in a straightforward manner. Apart from the drawback of restricting the generality to a certain class of ligament shapes included in the assumed model function, this has the advantage of dealing with the ligament as a whole dataset. Optimization strategies require parameter bounds, parameter start values of the model function. They are furthermore computationally demanding, because the parameter identification must be carried out individually and independently for each ligament. Therefore, we focus on the development of an artificial neural network (ANN).

For details on the ANNs, we refer to Huber (2018) and the literature cited there. For the training of the ANN, the 16 symmetric ligament geometries are used, defined in Table 1. Pattern files are written in the following style: The input vector X consists according to Equation (2) of the radii computed for all elements along the ligament skeleton, normalized by their average, in the form r(.)(xi)/〈r(.)〉, where (.) can be set to any of the three algorithms, namely (Th), (SE), or (EDT). This set of data represents the shape of the ligament from one end to the other end as measured by the corresponding algorithm.

As one further input value, the normalized position 2x_i/l is given, for which the correction factor shall be determined. The output vector Y consists according to Equation (3) of just one value, which is the correct radius divided by the radius determined from the algorithm r(ref)(xi)/r(.)(xi) at the position x_i. Because an ANN represents a continuous approximation of the presented data, it is very difficult to predict the steps contained in r^(ref) as shown in Figure 4. Therefore, the prediction of the output is limited to the positions within the triple points. Per ligament, 14 patterns are created, which are related to the 14 element radii for which the correct radius needs to be computed. Each ANN consists of four layers with 21 neurons at the input layer, 15 and 10 neurons in the two hidden layers, and 1 neuron for the output layer and is trained for 10,000 epochs. The resulting mean squared training and validation errors are MSET(Th)=2.37·10-5 and MSEV(Th)=1.87·10-4; MSET(SE)=8.82·10-6 and MSEV(SE)=8.44·10-5; MSET(EDT)=1.73·10-5 and MSEV(EDT)=1.89·10-4; respectively.

Indicated by the very low training and validation error, the reconstruction of the correct ligament shape seems to be a simple task for the ANN. The ANNs are able to determine the original ligament radius within one voxel accuracy, independent by which algorithm the input data are provided.

To validate the generalization capability of the trained ANNs, three new validation geometries are generated within the range of the existing 16 geometries. They are defined with r_end = 0.1375 and r_mid/r_end = [0.625, 0.875, 1.125]. The geometries of the three validation examples are shown in Figure 13. The degree of the remaining deviations is illustrated by 1 voxel- (±0.5 v) and 2 voxel-wide (±1 v) bands. The radii determined along the ligaments is within or very close to the 2-voxel wide band range for all three validation geometries and three ANN types. This corresponds to plus-minus one voxel, which is the limit for the accuracy defined by the voxel resolution. Only for the corrected Thickness data of strongly concave ligament in Figure 13A, the determined values are far outside the 2-voxel wide band.

If the Thickness data are pre-processed with the Smallest Ellipse algorithm, the accuracy improves significantly also for this difficult case. The ANN is now able to achieve accuracies, which are within the theoretical resolution limit of the voxel discretization and comparable to the results of the EDT (Figures 13A–C). This results from the capability of the ANN to memorize the relationship between ligament shapes and their corrections as whole and smoothly interpolate this relationship for untrained geometries. Due to this, the ANN approach has superior performance compared to the local Smallest Ellipse approach, discussed in the previous section, which cannot fully recover the information in the thinnest part of the ligament. The drawback is however that this method is so far limited to symmetric ligaments. For further evaluations of actual tomography data in the parameter space of rsym* and rasym*, as found by Richert and Huber (2018), the incorporation of a linear gradient according to Equation (1) is required. With this, asymmetric ligaments can be represented, as they occur with high probability in real NPG. Motivated by the promising results presented in this section, such an extension will be scope of future work.

As for the resulting mechanical behavior, very small deviations of maximum 10% are observed for the 16 trained geometries for both Young's modulus and yield strength. Also the three additional validation examples are well-predicted by the ANN with the very same accuracy, supposed the Thickness data are improved by the Smallest Ellipse approach before feeding the data to the ANN. It is remarkable that, despite some remaining error in the geometry reconstruction, resulting errors in the mechanical properties are negligible. The reason for this is that the ANN in average determines the ligament shape correctly with perhaps some small over- and under-predictions in different regions of the ligament. In contrast to that, the Thickness algorithm and the reconstruction via the Smallest Ellipse approach systematically overestimate the geometry and therefore the mechanical properties are biased to higher values.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.