Skeletonization, Geometrical Analysis, and Finite Element Modeling of Nanoporous Gold Based on 3D Tomography Data

: Various modeling approaches simplify and parametrize the complex network structure of nanoporous gold (NPG) for studying the structure–property relationship based on artiﬁcially generated structures. This paper presents a computational efﬁcient and versatile ﬁnite element method (FEM) beam model that is based on skeletonization and diameter information derived from the original 3D focused ion beam-scanning electron microscope (FIB-SEM) tomography data of NPG. The geometrical skeleton network is thoroughly examined for a better understanding of the NPG structure. A skeleton FEM beam model is derived that can predict the macroscopic mechanical behavior of the material. Comparisons between the mechanical response of this skeleton beam model and a solid FEM model are conducted. Results showed that the biggest-sphere diameter algorithm implemented in the open-source software FIJI, commonly used for geometrical analysis of microstructural data, overestimates the diameter of the curved NPG ligaments. The larger diameters lead to a signiﬁcant overestimation of macroscopic stiffness and strength by the skeleton FEM beam model. For a parabolic shaped ligament with only 20% variation in its diameter, a factor of more than two was found in stiffness. It is concluded that improved algorithms for image processing are needed that provide accurate diameter information along the ligament axis.


Introduction
Nanoporous gold (NPG) made by dealloying allows producing macroscopic millimeter-or even centimeter-sized porous bodies with a solid fraction around 30% [1][2][3][4]. The material exhibits a uniform, bicontinuous network of nanoscale pores and solid gold ligaments. The ligament diameter can be adjusted by annealing from a few nanometers to the micrometer scale [3,[5][6][7]. The mechanical behavior of NPG is a critical aspect of its functionality. Besides the strong and well-investigated size effect (see [2] and literature cited there), further work revealed an anomalously high elastic compliance, early yielding, and substantial accumulation of defects within the nanosized ligaments [6,[8][9][10][11][12].
To obtain deeper insights about the mechanical behavior and the underlying mechanisms under deformation, different modeling approaches are applied. Besides very fundamental calculations such as density functional theory (DFT) [13] and molecular dynamics simulations (MD) [11,12,14], there are several simplified finite element method (FEM) models [6,[13][14][15][16][17][18] and even analytical models [19][20][21][22][23][24] for predicting mechanical properties. In many modeling approaches the complex NPG Before the skeleton data were translated into the skeleton FEM beam model, the geometrical data were thoroughly analyzed. This helps with getting a sense of the skeleton and diameter data of the complex NPG structure generated by FIJI. The distribution of the diameter data as analyzed in detail by Hu et al. [7,15], gives a first impression about the overall structural characteristics. However, to create a realistic FEM beam model based on the skeleton and diameter information, the non-local information-i.e., the diameter evolution along the ligament-is also important. The object-oriented data processing with additionally programmed object functions makes this type of analysis possible. The skeleton structure including the dangling branches and the calculated diameters, both with special attention to boundary regions, are of particular interest. Furthermore, the whole set of connected ligaments will be analyzed with respect to the shape distribution. Such work was done previously by Pia et al. [28] based on 2D SEM images by visual fitting parabolic ligament shapes. Here, the 3D dataset provides full information and ligament shapes are not limited to a prior chosen approximating function.

Skeleton Branch Type Occurences and Length Distributions
In this section, the branches of the skeleton computed by FIJI are analyzed. Figure 1a shows a cutout from such a tomography 3D reconstruction together with the derived skeleton line and corresponding junction-and end-voxels, visualized as spheres colored in blue and red, respectively. It is noted that the skeleton branches span the whole connecting length to the node centers, because ligaments and nodes are reduced in their extension towards the center voxel. The solid ligaments would be reasonably shorter with a continuous transition into the nodal masses attached at both ends. This difference is important with regard to the expected mechanical response of the two modelling approaches, as it will be discussed later in Section 5.2.
The "AnalyzeSkeleton" plugin in FIJI allows to compute the distribution of branch lengths [36]. However, there is no differentiation between interconnected branches (branch connects two junctions), and dangling branches (branch connects junction with end). Looking at the dangling branches in Figure 1a, many seem to represent pinched-off ligaments formed during the annealing process, as discussed by Hu et al. [7,15]. Furthermore, some very short dangling branches (junctions and ends almost overlap) are present, which do not represent significant features of the structure. Additionally, some of the dangling branches end close to the RVE boundary as it is also schematically illustrated in Figure 1b. These artificial dangling branches (green) represent ligaments, which are cut by the RVE boundary and should be treated as interconnected branches. These should be included in the definition of the boundary conditions of the finite element beam model, while the real dangling branches should be free to move (red).
Metals 2017, 7, x FOR PEER REVIEW 4 of 20 illustrated in Figure 1b. These artificial dangling branches (green) represent ligaments, which are cut by the RVE boundary and should be treated as interconnected branches. These should be included in the definition of the boundary conditions of the finite element beam model, while the real dangling branches should be free to move (red).
(a) (b) Figure 1. (a) 3D cutout of the NPG tomography reconstruction [7,15]; skeleton-voxels are shown in black in original size; blue spheres mark the junction-voxels, red spheres mark the free end-voxels, and green spheres the end-voxels at the boundary; (b) Schematic 2D drawing of an RVE: The green circles illustrate the ligament ends that are due to the cuts at the RVE boundary. The gray shaped area represents the boundary region where these cut boundary ligaments were found. Due to the underlying method of surface thinning, the cut-ends do not always exactly lie at the surface, but further inside the RVE structure, as exemplified close to the right boundary.
For clarification of the number of occurrences of these boundary-dangling branches, a statistical analysis was conducted for all end-voxels. Therefore, the number of all end-voxels was counted depending on the distance from the nearest boundary plane, schematically shown as grey area in Figure 1b. A constant distribution of end-voxels was observed in the inner volume, but a four-fold number of end-voxels is detected in the boundary regions having a width of 3% of the cubic-length for all six boundary planes. This boundary layer width corresponds to 180 nm, which is very close to the average ligament radius of 200 nm. This suggests that the number of actual dangling branches (from local ligament thickenings and from pinched-off ligaments) is constant over the whole RVE volume. The additional end-voxels in the boundary region can be attributed to cut boundary ligaments, which can be assumed to be interconnected to the structure outside of the RVE. It should be noted that the counts are evenly spread in this boundary region, i.e. most cut-branches are not exactly located on the RVE boundary. This is due to the method of surface thinning that is used for skeletonization. Thus, end-voxels within a distance of the average ligament radius from the boundary region need to be included in the boundary conditions of the model.
In Figure 2, the stacked histogram illustrates the distribution of the branch length and occurring frequencies for the three types of branches. The blue bins represent the frequency of occurrence of fully-interconnected branches. The red bins add the frequency of free-dangling branches and the green ones the boundary-dangling ligaments lying in the narrow boundary region. The total sum represents the histogram for all skeleton-branches. 912 of the total 1284 skeletonbranches (71%) are fully-interconnected within the RVE, 208 are free-dangling branches (16%), and 164 are boundary-dangling branches (13%). This furthermore shows that the boundary-dangling branches need to be constrained correctly to mimic realistic boundary conditions. The mean freedangling branch length is calculated as 329 nm, which is again slightly smaller than the mean ligament diameter of 400 nm of the structure. The mean length of the boundary-dangling branches is  [7,15]; skeleton-voxels are shown in black in original size; blue spheres mark the junction-voxels, red spheres mark the free end-voxels, and green spheres the end-voxels at the boundary; (b) Schematic 2D drawing of an RVE: The green circles illustrate the ligament ends that are due to the cuts at the RVE boundary. The gray shaped area represents the boundary region where these cut boundary ligaments were found. Due to the underlying method of surface thinning, the cut-ends do not always exactly lie at the surface, but further inside the RVE structure, as exemplified close to the right boundary.
For clarification of the number of occurrences of these boundary-dangling branches, a statistical analysis was conducted for all end-voxels. Therefore, the number of all end-voxels was counted depending on the distance from the nearest boundary plane, schematically shown as grey area in Figure 1b. A constant distribution of end-voxels was observed in the inner volume, but a four-fold number of end-voxels is detected in the boundary regions having a width of 3% of the cubic-length for all six boundary planes. This boundary layer width corresponds to 180 nm, which is very close to the average ligament radius of 200 nm. This suggests that the number of actual dangling branches (from local ligament thickenings and from pinched-off ligaments) is constant over the whole RVE volume. The additional end-voxels in the boundary region can be attributed to cut boundary ligaments, which can be assumed to be interconnected to the structure outside of the RVE. It should be noted that the counts are evenly spread in this boundary region, i.e., most cut-branches are not exactly located on the RVE boundary. This is due to the method of surface thinning that is used for skeletonization. Thus, end-voxels within a distance of the average ligament radius from the boundary region need to be included in the boundary conditions of the model.
In Figure 2, the stacked histogram illustrates the distribution of the branch length l and occurring frequencies for the three types of branches. The blue bins represent the frequency of occurrence of fully-interconnected branches. The red bins add the frequency of free-dangling branches and the green ones the boundary-dangling ligaments lying in the narrow boundary region. The total sum represents the histogram for all skeleton-branches. 912 of the total 1284 skeleton-branches (71%) are fully-interconnected within the RVE, 208 are free-dangling branches (16%), and 164 are boundary-dangling branches (13%). This furthermore shows that the boundary-dangling branches need to be constrained correctly to mimic realistic boundary conditions. The mean free-dangling branch length is calculated as 329 nm, which is again slightly smaller than the mean ligament diameter of 400 nm of the structure. The mean length of the boundary-dangling branches is 539 nm and the mean length of the fully-interconnected branches is 692 nm. An asymmetric distribution is visible in the histogram for all branch types.  A large number of 105 free-dangling branches is visible at shortest branch length (first red bin), which corresponds to half of all free-dangling branches. These short free-dangling branches are a general uncertainty of skeletonization algorithms and result from the sensitivity of the algorithm to noise at the surface or in general to the ligament shape. A local thickening of a ligament typically leads to an extra junction with a very short dangling branch, which is fully enclosed by the solid ligament. The longer free-dangling branches are assumed to be actual pinched-off ligaments formed during the annealing process of the sample [7,15]. For a skeleton pruning, a criterion for differentiation between unwanted (structurally not meaningful) and pinched-off ligaments would be needed. However, the dangling branches do not influence the mechanical response of the simulation, because they do not carry load. Thus, the development of such a criterion is not relevant for the further course of this work.

Skeleton-Based Diameter Distribution and Regional Occurance
In addition to the skeleton data, the diameter information estimated by the "Thickness" algorithm is needed for the skeleton FEM beam model. Hu et al. thoroughly investigated the diameter results of the volume-based diameter data calculated by FIJI with respect to the self-similarity of the NPG structure [7,15]. "Volume-based" implies that the diameter analysis by FIJI is done for all voxels within the solid phase. Our model is based on the skeleton and thus just takes the diameter information attributed to the skeleton-voxels into consideration. To assure comparability, the mean ligament diameter, as well as the diameter distribution, are determined for the volume-based and skeleton-based approach, respectively. For the volume-based approach, a mean diameter 〈 〉 of 414 nm is calculated, while a slightly smaller (3%) mean ligament diameter 〈 〉 of 401 nm results for the skeleton-based approach. Furthermore, Hu et al. showed for several RVE sizes that the volume-based diameter distribution represents a Gaussian distribution [7,15]. The skeleton-based diameter distribution can also be well represented by the Gaussian distribution, which is more flat and slightly shifted to smaller diameters though.
The detected amount of smaller diameters that appear in the skeleton-based diameter distribution motivated a refined analysis in dependence of the position within the RVE. The question arises if these skeleton-voxels are located in any specific region of the RVE structure. To this end, the whole RVE was divided into four regions with constant distance from all boundaries, all having the same thickness of 30 voxels. The outermost region was again subdivided considering the narrow A large number of 105 free-dangling branches is visible at shortest branch length (first red bin), which corresponds to half of all free-dangling branches. These short free-dangling branches are a general uncertainty of skeletonization algorithms and result from the sensitivity of the algorithm to noise at the surface or in general to the ligament shape. A local thickening of a ligament typically leads to an extra junction with a very short dangling branch, which is fully enclosed by the solid ligament. The longer free-dangling branches are assumed to be actual pinched-off ligaments formed during the annealing process of the sample [7,15]. For a skeleton pruning, a criterion for differentiation between unwanted (structurally not meaningful) and pinched-off ligaments would be needed. However, the dangling branches do not influence the mechanical response of the simulation, because they do not carry load. Thus, the development of such a criterion is not relevant for the further course of this work.

Skeleton-Based Diameter Distribution and Regional Occurance
In addition to the skeleton data, the diameter information estimated by the "Thickness" algorithm is needed for the skeleton FEM beam model. Hu et al. thoroughly investigated the diameter results of the volume-based diameter data calculated by FIJI with respect to the self-similarity of the NPG structure [7,15]. "Volume-based" implies that the diameter analysis by FIJI is done for all voxels within the solid phase. Our model is based on the skeleton and thus just takes the diameter information attributed to the skeleton-voxels into consideration. To assure comparability, the mean ligament diameter, as well as the diameter distribution, are determined for the volume-based and skeleton-based approach, respectively. For the volume-based approach, a mean diameter d of 414 nm is calculated, while a slightly smaller (3%) mean ligament diameter d of 401 nm results for the skeleton-based approach. Furthermore, Hu et al. showed for several RVE sizes that the volume-based diameter distribution represents a Gaussian distribution [7,15]. The skeleton-based diameter distribution can also be well represented by the Gaussian distribution, which is more flat and slightly shifted to smaller diameters though. The detected amount of smaller diameters that appear in the skeleton-based diameter distribution motivated a refined analysis in dependence of the position within the RVE. The question arises if these skeleton-voxels are located in any specific region of the RVE structure. To this end, the whole RVE was divided into four regions with constant distance from all boundaries, all having the same thickness of 30 voxels. The outermost region was again subdivided considering the narrow boundary region with a width of 3% of the cubic length (seven voxels), found in the previous Section 3.1. The inner region with a voxel-distance of 90-120 to any boundary has the smallest volume, as schematically illustrated by the cross-sectional cut through the RVE in Figure 3a. For comparability, Gaussian curves were fit to the diameter distributions of the selected regions and normalized to an area of one, plotted in Figure 3b.  For comparison, a curve of the diameter distribution of all skeleton voxels in the whole RVE is plotted as a dashed black line. It is clearly visible that only the distribution for the outermost narrow 0-7 region shows a left shift to smaller ligament diameters and is more flat. This clearly shows the existence of non-representative small diameters, which exists only in the boundary region. The distribution for the whole RVE just shows a slight influence by this boundary region.
The smaller diameters located at the RVE boundary are caused by the cutting of the ligaments at the boundary, as it is exemplarily visible from the RVE in Figure 1a for the ligaments at the right boundary. Inclined ligaments are falsely interpreted as increasingly thinned as they approach the RVE boundary, and are receiving artificially low diameter information. This is important for the boundary conditions applied to the skeleton FEM beam model and is further discussed in Section 4.2.

Ligament Shape Distribution
For the shape analysis of the NPG structure, most publications differentiate between nodes (or nodal masses) and ligaments. They consider the nodes as spherical [6,18,23] or cubic [25,28] volume extending from the connecting junctions, and the ligaments as the volume between these nodes. In this work, the nodal volume corresponds to the biggest sphere containing the according junctionvoxel, as it is fit into the structure by the FIJI "Thickness" plugin [38]. For the following translation into a finite element model (see Section 4), beam elements will be created along the skeleton-line voxels considering the attributed diameters. Thus, the shape along the whole branch length matters for the mechanical behavior and is analyzed as part of the ligaments in the following. Furthermore, to ensure a clean analysis, we consider only the fully interconnected branches and of those only branches that are equal or longer than the mean value of both node (junction) diameters. This selection excludes branches where the nodal volumes are too close and overlap. 793 branches of the overall 1284 skeleton branches remain.
Concerning the parametrization of the contour of the ligaments, several papers state that a For comparison, a curve of the diameter distribution of all skeleton voxels in the whole RVE is plotted as a dashed black line. It is clearly visible that only the distribution for the outermost narrow 0-7 region shows a left shift to smaller ligament diameters and is more flat. This clearly shows the existence of non-representative small diameters, which exists only in the boundary region. The distribution for the whole RVE just shows a slight influence by this boundary region.
The smaller diameters located at the RVE boundary are caused by the cutting of the ligaments at the boundary, as it is exemplarily visible from the RVE in Figure 1a for the ligaments at the right boundary. Inclined ligaments are falsely interpreted as increasingly thinned as they approach the RVE boundary, and are receiving artificially low diameter information. This is important for the boundary conditions applied to the skeleton FEM beam model and is further discussed in Section 4.2.

Ligament Shape Distribution
For the shape analysis of the NPG structure, most publications differentiate between nodes (or nodal masses) and ligaments. They consider the nodes as spherical [6,18,23] or cubic [25,28] volume extending from the connecting junctions, and the ligaments as the volume between these nodes. In this work, the nodal volume corresponds to the biggest sphere containing the according junction-voxel, as it is fit into the structure by the FIJI "Thickness" plugin [38]. For the following translation into a finite element model (see Section 4), beam elements will be created along the skeleton-line voxels considering the attributed diameters. Thus, the shape along the whole branch length matters for the mechanical behavior and is analyzed as part of the ligaments in the following. Furthermore, to ensure a clean analysis, we consider only the fully interconnected branches and of those only branches that are equal or longer than the mean value of both node (junction) diameters. This selection excludes branches where the nodal volumes are too close and overlap. 793 branches of the overall 1284 skeleton branches remain.
Concerning the parametrization of the contour of the ligaments, several papers state that a minimum diameter can be observed roughly in the middle, with increasing diameter approaching the nodes [28,41]. Pia et al. [28] analyzed SEM images by fitting a parabolic function of the form to the ligament shape visible in 2D images.
The shape factor b gives information about the concaveness (b > 0) or convexness (b < 0) of the ligaments. Pia et al. though just consider concave ligaments. This idea furthermore, as also in other literature, assumes the shape of the ligaments to be symmetrical. Towards a more detailed investigation, the object-oriented data structure established in this work conveniently allows for a more general statistical analysis of the ligament shape. We additionally consider ligaments that have an asymmetric contour, besides the known symmetric concave, convex, and cylindrical shapes. For the analysis, we divide the ligament into three segments of equal length: two end and one middle segment. From the position and diameter information stored in the skeleton voxel objects, the average radii are calculated as r 1 and r 2 for both end segments, respectively, and r m for the middle segment. The results for the three segments are organized such that r 1 ≤ r 2 applies for all ligaments. A ligament is cylindrical when the radius is constant (r 1 = r 2 = r m ), it is convex for r 1 , r 2 < r m or concave for r 1 , r 2 > r m . As a particular characteristic among the asymmetric shapes, the s-shape contours follow the condition: For an overview of all ligament shapes, two independent shape parameters are defined, which characterize the symmetric convexity/concavity and the asymmetry of a ligament.
These two independent shape parameters are intended to give a first qualitative and statistically robust overview of the distribution of ligament shapes within the whole RVE. The statistical distribution of r * sym and r * asym are shown in the point-plot in Figure 4. The left region includes concave ligaments, the middle the s-shaped ligaments (green points) and the right region the convex shapes. We account all ligaments as cylindrical, where r * sym or the asymmetry factor r * asym is up to ±10%. These are found in the red box around the point r * sym = 1 and r * asym = 0 with a size ±0.1. The green lines at ±63.4 • angle indicate the transition into the s-shape. At these lines, the average radius of either end-segment r 1 or r 2 equals the radius of the middle segment, r m . The inserts in Figure 4 give an impression of the ligament shapes in dependence of the position in the diagram. The total of all data points within the corresponding regions is given in Table 1, delivering a quantitative overview of the probability of the different shapes.
It turns out that 200 ligaments are concave (r 1 , r 2 > r m ) in comparison to just 31 convex (r 1 , r 2 < r m ) ligaments. 22 of the initially 53 detected convex ligaments are counted as cylindrical according to the 10% tolerance. Cylindrical ligaments appear with about the same frequency as the concave ligaments. However, this number strongly depends on the chosen cut-off for the definition of the cylindrical shape. Furthermore, a surprisingly high number of 370 s-shaped ligaments is observed.
High asymmetry values are observed for such ligament shapes, representing ligaments that are connecting a very thin and a very thick node.
For obtaining a better overview on the distribution of the shape characteristics, Figure 5a presents a histogram revealing the distribution of the symmetric properties. Considering only the convex, concave, and cylindrical ligaments included in the ±10% tolerance, the majority of ligaments are found in a range r * sym from 0.7 to 1.1 with a maximum close to 1.0. The histogram in Figure 5b shows the asymmetry r * asym . Up to an asymmetry ratio of 0.1, all ligaments are convex, concave or, due to the chosen tolerance, cylindrical. From 0.1 up to 0.55 the frequency of total ligaments drops, with increasing proportion of s-shaped ligaments. At r * asym = 0.1 a significant fraction of at least 50% originates from s-shaped ligaments. Only a few ligaments show an asymmetry higher than 0.7.  . Ligament shape analysis: The ligaments were divided into three equal length segments with the average radii and for both end segments, and for the middle segment. Point-plot of cylindrical, concave, and convex (CCC) (blue points) and s-shape (green) contoured ligaments in dependence of the two independent shape parameters * and * , defining the degree of concavity ( * < 1), or convexity ( * > 1), and asymmetry, respectively. Cylindrical ligaments are defined with a tolerance of ±0  . Ligament shape analysis: The ligaments were divided into three equal length segments with the average radii r 1 and r 2 for both end segments, and r m for the middle segment. Point-plot of cylindrical, concave, and convex (CCC) (blue points) and s-shape (green) contoured ligaments in dependence of the two independent shape parameters r * sym and r * asym , defining the degree of concavity (r * sym < 1), or convexity (r * sym > 1), and asymmetry, respectively. Cylindrical ligaments are defined with a tolerance of ±0.1, represented by the red box around the point r * sym = 1 and r * asym = 0. The inserted schematic 3D drawings show exemplary ligament contours at the points marked by the red circles. . Ligament shape analysis: The ligaments were divided into three equal length segments with the average radii and for both end segments, and for the middle segment. Point-plot of cylindrical, concave, and convex (CCC) (blue points) and s-shape (green) contoured ligaments in dependence of the two independent shape parameters * and * , defining the degree of concavity ( * < 1), or convexity ( * > 1), and asymmetry, respectively. Cylindrical ligaments are defined with a tolerance of ±0.1, represented by the red box around the point * = 1 and * = 0. The inserted schematic 3D drawings show exemplary ligament contours at the points marked by the red circles.  with concave, convex, or cylindrical shape (CCC) and s-shape in dependence of the parameters: (a) r * sym characterizing the degree of concavity (r * sym < 1) or convexity (r * sym > 1); (b) r * asym characterizing the degree of asymmetry. This analysis highlights the diversity of ligament shapes that are found in the NPG structure. Almost half of all analyzed ligaments have an asymmetric and even s-shaped contour. Thus, models limited to symmetrical ligament geometries-e.g., parabolic shapes-are insufficient. In conclusion, this analysis underlines the importance of considering the actual ligament shapes from 3D tomography data for the prediction of a realistic mechanical response of NPG. This is content of the following Section 4, which aims on translating the skeleton and diameter information into a computationally efficient finite element beam model.

Skeleton FEM Beam Model
The FEM model to be developed in this section should incorporate the real microstructure but, at the same time, it should allow also for a high computational efficiency. To this end, the discretization of the structure will make use of beam elements available in ABAQUS [40]. As proposed by Huber et al. [6,[16][17][18], cylindrical B31 beam elements are used, which allow transverse shear strain considering Timoshenko beam theory.
The novelty in comparison to the artificially generated beam model of Huber et al. is, that the skeleton FEM beam model emulates the actual ligament path and shape of the NPG sample, based on the skeleton and diameter information, respectively, extracted from the NPG tomography. As for the translation of the tomography data into the beam model, each ligament is divided into several beam elements along the skeleton line. The number of beam elements per ligament depends on the chosen discretization, which is discussed in Section 4.1. Each cylindrical beam element receives the average diameter of the associated voxels.

Discretization of Branches
To study the influence of different possible discretization approaches, the same cutout of the RVE is illustrated in Figure 6b-d for each skeleton FEM beam model, and discussed in the following. The cutout of the corresponding solid FEM model created from the same tomography data reconstruction of Hu et al. [7,15] is represented in Figure 6a. Besides the different number of beam elements per branch, the different beam diameters are clearly visible in Figure 6b-d, representing the ligament contour.
For a voxelized skeleton, the two extreme discretization approaches are: (i) using just one beam per branch, resulting in a straight connection between two junctions, or (ii) representing the connection between two neighbor voxels by an individual beam element (1 V/E). The first approach (i) was originally used by Hu et al. [7,15] to investigate the idea of load-bearing ring structures. It is clear that straight connections, as shown in Figure 6b, will lead to an overestimation in stiffness and strength. The curved geometry of the initial structure is not represented. This gets obvious for the ring that is present in the initial structure, which is quasi narrowed down to beams along its center axis with the discretization (i). Furthermore, the ligaments that actually first meander out of the cutout volume and back in are as straight connections already displayed within the cutout volume. In contrast to this discretization, the second approach (ii) would lead to a zigzag-like structure with strongly changing relative axis orientation of, e.g., 35 • , 45 • , or 90 • between the beam elements, as defined by the next neighbor in the grid of cubic voxels, see Figure 6d. In this case, the mechanical response is expected to be too compliant and too soft. A reasonable discretization should lie somewhere in between.
Due to the broad distribution in branch length, as discussed in Section 3.1 (Figure 2) it is not sufficient to divide every branch into the same number of beam elements. For short branches, this would lead to artificial FE nodes that lie between the voxels; for very long branches, the discretization might be insufficient for resolving the path of the ligament. Therefore, the implemented discretization approach aims for an average beam element length of five voxels. All branches shorter than five voxels will be represented by one single beam element. Longer branches will be divided into beam elements that each span five voxels, if the number of voxels is a multiple of five. For all other cases, the element length is varied between four or six voxels, depending on the branch length. This leads to a reasonable representation of branches in the skeleton FEM beam model, referred to as 5 V/E in the following, see Figure 6c. Further visual comparisons of planar thin boundary sections of the two modeling approaches were carried out, showing the same distinctive features and thus validating the geometrical matching.
With the chosen 5 V/E discretization approach, the average number of elements per fully-interconnected branch in the RVE is 7. In contrast to that, Jiao et al. found that approximately 10 elements per branch are necessary to predict the mechanical behavior of a ligament with sufficient accuracy [17]. Further refinement is, however, limited in our model because of the number of voxels per ligament. The 5 V/E discretization therefore represents a compromise regarding a reasonable smoothening based on the typical number of voxels per ligament and desired computational accuracy of the model. approach aims for an average beam element length of five voxels. All branches shorter than five voxels will be represented by one single beam element. Longer branches will be divided into beam elements that each span five voxels, if the number of voxels is a multiple of five. For all other cases, the element length is varied between four or six voxels, depending on the branch length. This leads to a reasonable representation of branches in the skeleton FEM beam model, referred to as 5 V/E in the following, see Figure 6c. Further visual comparisons of planar thin boundary sections of the two modeling approaches were carried out, showing the same distinctive features and thus validating the geometrical matching.
With the chosen 5 V/E discretization approach, the average number of elements per fullyinterconnected branch in the RVE is 7. In contrast to that, Jiao et al. found that approximately 10 elements per branch are necessary to predict the mechanical behavior of a ligament with sufficient accuracy [17]. Further refinement is, however, limited in our model because of the number of voxels per ligament. The 5 V/E discretization therefore represents a compromise regarding a reasonable smoothening based on the typical number of voxels per ligament and desired computational accuracy of the model.

Material Properties and Boundary Conditions
All simulations were performed with the following material parameters, being the elastic modulus of gold = 81 GPa, Poisson's ratio = 0.42, yield stress = 700 MPa and work hardening rate = 1000 MPa . These values are taken from Hu et al. [7,15] and allow later comparison of the predicted macroscopic properties from both simulations, i.e., skeleton vs. solid model.
In the solid FEM model of Hu et al. [7,15], symmetric boundary conditions were applied to all FE nodes on the three planes = 0, = 0, and = 0. These FE nodes belong to the cut-surfaces of the NPG ligaments, which would continue beyond the boundary in the larger specimen. All other FE nodes on the remaining side faces were free to move. The load was applied as a homogeneous displacement of all FE nodes on the top side of the RVE = . For comparability, equivalent boundary conditions are applied for the skeleton FEM beam model. All end voxels in the 3% cubic length boundary regions are assumed to be cut boundary-ligaments, as discussed in Section 3.3. Furthermore, the "Thickness" algorithm produces considerably smaller diameter information for inclined ligaments in the boundary region. Thus, symmetric boundary conditions are applied to all FE nodes in this 3% wide boundary region at the three planes = 0, = 0, and = 0. This corresponds to only 10% of the total volume. For applying the load, all FE nodes selected in the top boundary region = are displaced by 8% macroscopic compression strain, which results in approximately 477 nm compression displacement.

Material Properties and Boundary Conditions
All simulations were performed with the following material parameters, being the elastic modulus of gold E Au = 81 GPa, Poisson's ratio ν Au = 0.42, yield stress σ y = 700 MPa and work hardening rate γ Au = 1000 MPa. These values are taken from Hu et al. [7,15] and allow later comparison of the predicted macroscopic properties from both simulations, i.e., skeleton vs. solid model.
In the solid FEM model of Hu et al. [7,15], symmetric boundary conditions were applied to all FE nodes on the three planes x = 0, y = 0, and z = 0. These FE nodes belong to the cut-surfaces of the NPG ligaments, which would continue beyond the boundary in the larger specimen. All other FE nodes on the remaining side faces were free to move. The load was applied as a homogeneous displacement of all FE nodes on the top side of the RVE z = z max . For comparability, equivalent boundary conditions are applied for the skeleton FEM beam model. All end voxels in the 3% cubic length boundary regions are assumed to be cut boundary-ligaments, as discussed in Section 3.3. Furthermore, the "Thickness" algorithm produces considerably smaller diameter information for inclined ligaments in the boundary region. Thus, symmetric boundary conditions are applied to all FE nodes in this 3% wide boundary region at the three planes x = 0, y = 0, and z = 0. This corresponds to only 10% of the total volume. For applying the load, all FE nodes selected in the top boundary region z = z max are displaced by 8% macroscopic compression strain, which results in approximately 477 nm compression displacement.

Mechanical Response of the Skeleton FEM Beam Model
The mechanical responses predicted by the new skeleton FEM beam model were analyzed in comparison to the FEM simulations based on the solid FEM model consisting of 10-node tetrahedral elements (C3D10) provided by Hu [7]. In the following, only the compressive loading in z-direction is discussed. The compression of the RVE in x-and y-direction showed similar anisotropy, as Hu described in [7]. The major advantage of the skeleton FEM beam model becomes obvious by comparing the computation times of 0.022 h for the skeleton beam model and 1.38 h for the solid model for elastic-plastic material behavior and same total strain.

Simulation Results for RVE
The skeleton FEM beam model with the 5 V/E discretization and the corresponding solid FEM model are illustrated in Figure 7a-b, respectively. The color bar represents the von Mises stress ranging from 0 MPa to a user-defined maximum of 1000 MPa. With this stress color coding, locally higher stress concentrations are visible in the solid FEM model. The skeleton FEM beam model does not show as many high stress concentrations but a more evenly distributed stress field. Nevertheless, the overall stress distribution appears to be comparable. The mechanical responses predicted by the new skeleton FEM beam model were analyzed in comparison to the FEM simulations based on the solid FEM model consisting of 10-node tetrahedral elements (C3D10) provided by Hu [7]. In the following, only the compressive loading in z-direction is discussed. The compression of the RVE in x-and y-direction showed similar anisotropy, as Hu described in [7]. The major advantage of the skeleton FEM beam model becomes obvious by comparing the computation times of 0.022 h for the skeleton beam model and 1.38 h for the solid model for elastic-plastic material behavior and same total strain.

Simulation Results for RVE
The skeleton FEM beam model with the 5 V/E discretization and the corresponding solid FEM model are illustrated in Figure 7a-b, respectively. The color bar represents the von Mises stress ranging from 0 MPa to a user-defined maximum of 1000 MPa. With this stress color coding, locally higher stress concentrations are visible in the solid FEM model. The skeleton FEM beam model does not show as many high stress concentrations but a more evenly distributed stress field. Nevertheless, the overall stress distribution appears to be comparable. . Stress distribution in the RVE produced from the same NPG structure for the two simulation approaches, both compressed to 8% macroscopic strain, the legend shows the von Mises stress with the color maximum red set to a stress of 1000 MPa: (a) solid FEM model [7], (b) skeleton FEM beam model.

Macroscopic Mechanical Response
The resulting macroscopic stress-strain behavior was calculated for both modeling approaches and the stress-strain curves are plotted in Figure 8. For the skeleton beam approach, also the two extreme discretization approaches were calculated to give an impression on the effect of the discretization. The skeleton FEM beam model with approximately five voxels per beam element (5 V/E) is additionally tested with quadratic B32 beam elements as alternative to a finer discretization with 10 elements. The elastic modulus was determined from the slope of the initial linear part of the loading segment. It is difficult to define a precise yield point of NPG due to the continuous transition from elastic to plastic behavior that happens already at very small loads. Here, the offset plastic strength was defined at 0.2% plastic strain, as in [7,15]. The calculated results for all models are collected in Table 2.
The skeleton beam model with straight beam connections shows the highest Young's modulus and yield strength. This proves the expectation that the straight connections over-stiffen the network, in comparison to a network with curved connection paths. This is consistent to the finding of Huber et al. [6] and Roschning et al. [16] that the amount of disorder going along with increasing beam curvedness has a strong impact on the mechanical response of such a beam network. This furthermore goes along with the observation that the 1 V/E discretization shows highest compliance and the flattest stress-strain curve. This is due to the large angles between the individual beam elements, leading to a less stiff and less strong response. Furthermore, the simulation with this model Figure 7. Stress distribution in the RVE produced from the same NPG structure for the two simulation approaches, both compressed to 8% macroscopic strain, the legend shows the von Mises stress with the color maximum red set to a stress of 1000 MPa: (a) solid FEM model [7], (b) skeleton FEM beam model.

Macroscopic Mechanical Response
The resulting macroscopic stress-strain behavior was calculated for both modeling approaches and the stress-strain curves are plotted in Figure 8. For the skeleton beam approach, also the two extreme discretization approaches were calculated to give an impression on the effect of the discretization. The skeleton FEM beam model with approximately five voxels per beam element (5 V/E) is additionally tested with quadratic B32 beam elements as alternative to a finer discretization with 10 elements. The elastic modulus was determined from the slope of the initial linear part of the loading segment. It is difficult to define a precise yield point of NPG due to the continuous transition from elastic to plastic behavior that happens already at very small loads. Here, the offset plastic strength was defined at 0.2% plastic strain, as in [7,15]. The calculated results for all models are collected in Table 2.
The skeleton beam model with straight beam connections shows the highest Young's modulus and yield strength. This proves the expectation that the straight connections over-stiffen the network, in comparison to a network with curved connection paths. This is consistent to the finding of Huber et al. [6] and Roschning et al. [16] that the amount of disorder going along with increasing beam curvedness has a strong impact on the mechanical response of such a beam network. This furthermore goes along with the observation that the 1 V/E discretization shows highest compliance and the flattest stress-strain curve. This is due to the large angles between the individual beam elements, leading to a less stiff and less strong response. Furthermore, the simulation with this model was terminated at around 4% strain, which is probably due to occurring localizations and instabilities. The skeleton beam model with 5V/E, shows a stress-strain curve that lies between these two extreme cases, with an elastic modulus of about 480 MPa and a plastic strength of 8.0 MPa. Using the more accurate B32 beam elements, the difference to the result with B31 elements is negligible compared to the effect of the discretization. This indicates, that five elements per branch are sufficiently accurate compared to other modeling uncertainties.
In comparison to the yield strength of the solid model, all skeleton beam models show higher yield strengths by 40% (1 V/E), 105% (5 V/E), and 240% (straight connection). Also, the Young's modulus is higher by 30% (5 V/E) and 57% (straight connection). Only for 1 V/E, the Young's modulus coincides with that of the solid model. It is noted that in all cases the overestimation of the yield strength is higher than the overestimation of the elastic modulus.
Metals 2017, 7, x FOR PEER REVIEW 12 of 20 was terminated at around 4% strain, which is probably due to occurring localizations and instabilities. The skeleton beam model with 5V/E, shows a stress-strain curve that lies between these two extreme cases, with an elastic modulus of about 480 MPa and a plastic strength of 8.0 MPa. Using the more accurate B32 beam elements, the difference to the result with B31 elements is negligible compared to the effect of the discretization. This indicates, that five elements per branch are sufficiently accurate compared to other modeling uncertainties.
In comparison to the yield strength of the solid model, all skeleton beam models show higher yield strengths by 40% (1 V/E), 105% (5 V/E), and 240% (straight connection). Also, the Young's modulus is higher by 30% (5 V/E) and 57% (straight connection). Only for 1 V/E, the Young's modulus coincides with that of the solid model. It is noted that in all cases the overestimation of the yield strength is higher than the overestimation of the elastic modulus.  The results obtained from the skeleton FEM beam model are unexpected. Instead of an overprediction of stiffness and strength values lower than those from the solid model were expected, at least for the finer discretization 5 V/E and even more for 1 V/E. The reason for this expectation is the stiffening and strengthening effect of the nodal mass that is accounted for in the solid model and still missing in the skeleton beam model (see [6,16,18]).
In the beam model, the beam elements of different ligaments connect in a material point (vertex), and not in a massive and thus comparably stiff and strong node as present in the real NPG structure.  The results obtained from the skeleton FEM beam model are unexpected. Instead of an over-prediction of stiffness and strength values lower than those from the solid model were expected, at least for the finer discretization 5 V/E and even more for 1 V/E. The reason for this expectation is the stiffening and strengthening effect of the nodal mass that is accounted for in the solid model and still missing in the skeleton beam model (see [6,16,18]).
In the beam model, the beam elements of different ligaments connect in a material point (vertex), and not in a massive and thus comparably stiff and strong node as present in the real NPG structure. This nodal mass leads to a shortening of the actual ligament length identified as an important effect, mostly for the plastic behavior [6]. Jiao et al. conducted a study comparing a solid model with a beam model of a tetrahedron structure to calculate the necessary nodal correction in the beam model [18]. Their approach is to artificially strengthen the beam elements in the nodal area, mainly by increasing the diameter and extension of the elements making up the nodal region. It was concluded that, in this way, the Young's modulus and yield strength of a beam model with 30% solid fraction would be increased by around 120% and 80%, respectively, to match the values of the corresponding solid model [18]. Together with the results presented in Table 2, this adds up to an overestimation of almost a factor of two with the SBM 5 V/E B31 model.

Impact of "Thickness" Algorithm on Ligament Stiffness
From the previous section, it was concluded that macroscopic response of the skeleton beam model overestimates stiffness and strength, although the opposite was expected in comparison to a FEM solid model. This unexpected result suggests that there exists a common origin with regard to both properties. Besides the structure of the skeleton, which was carefully validated, see Figure 6, the ligament diameter is another critical parameter that has a strong effect on the mechanical response. The question arises, how well the "Thickness" algorithm [38,39] reproduces the real diameter in each point of the solid phase. The "Thickness" algorithm fits, for each voxel, the biggest possible sphere within the solid that contains this particular voxel. The sphere diameter is attributed to that individual voxel. This approach is often used for analyzing diameter distributions of fibers and bone structures [42], which are characterized by a nearly cylindrical shape. However, as can be seen from Table 1, 70% of the ligaments show strongly changing radii along the ligament axis by more than 10%. In contrast to the average ligament, which is indeed nearly cylindrical, the data for individual ligaments, see Figure 4, reveal that the radius variation along the ligament can be significant. The shape analysis revealed that a high percentage of ligaments are concave and possess symmetry factors up to r * sym = 0.5.
To study the output of the "Thickness" algorithm for a curved ligament, we exemplary investigate the radius results for a simple parabolic ligament with r * sym values in this observed range. Furthermore, we give an estimation about the effect on the predicted mechanical stiffness of a ligament with such curved contours. An exemplary concave ligament is shown in blue in the schematic drawing in Figure 9, parametrized according to Equation (1).
The ligament is assumed to transition into the nodal region at point Q, as it is tangentially touching the nodal spherical volume with the midpoint N. For the following analysis, the ligament length l is defined as the node-to-node-distance according to [6,16]. This is motivated by the skeleton FEM beam model, which connects two node center points using beam elements generated along the ligament axis.
For an arbitrarily chosen point x P on the ligament axis, the biggest sphere that contains this point has a radius R and midpoint M. It tangentially touches the boundary of the ligament at point T and is significantly shifted to a larger ligament radius than r(x P ). Following the "Thickness" algorithm, the radius R is now attributed as size information to the point x P . Repetition of this construction for all points of the ligament axis leads to the red curve r (x) entered in Figure 9, which corresponds to the identified shape of the original blue ligament with r(x). The comparison reveals that the radius of the ligament is systematically overestimated by an amount that increases with the curvature of the ligament. With transitioning into the nodal area in point Q, the identified radius takes the value of the spherical node, which results in a plateau until point N. factors up to = 0.5.
To study the output of the "Thickness" algorithm for a curved ligament, we exemplary investigate the radius results for a simple parabolic ligament with * values in this observed range. Furthermore, we give an estimation about the effect on the predicted mechanical stiffness of a ligament with such curved contours. An exemplary concave ligament is shown in blue in the schematic drawing in Figure 9, parametrized according to Equation (1). Figure 9. Schematic drawing of a parabolic ligament ( ) with length (blue curve), and construction of the identified shape ( ) according to the "Thickness" algorithm (red curve). The biggest sphere for the point has the midpoint . The ligament is connected to a spherical node with midpoint (grey shaded area). Figure 9. Schematic drawing of a parabolic ligament r(x) with length l (blue curve), and construction of the identified shape r (x) according to the "Thickness" algorithm (red curve). The biggest sphere for the point x P has the midpoint M. The ligament is connected to a spherical node with midpoint N (grey shaded area).
To calculate the actual quantitative radius overestimation and the resulting stiffness overestimation, a mathematical study was conducted for the parabolic shaped ligament. For the geometrical construction of the result of the "Thickness" algorithm for a given ligament geometry with parameters a and b, according to Equation (1), it is convenient to start from an arbitrary chosen point T = (x T , r T ), see Figure 9. It defines the position on the ligament axis x = x T , where the fitted sphere with radius R and center point M tangentially touches the ligament contour r(x). The line g(x) is constructed normal to the tangent in point T and cuts through M: Solving Equation (4) for g(x = x M ) = 0 yields Equation (5) from which the radius of the sphere R can be calculated using Equation (6). Finally, the corresponding position x P to which the radius R is assigned by the "Thickness" algorithm is determined from Equation (7): Using Equations (4)-(7) the shape of the ligament r (x P ) as determined by the "Thickness" algorithm can be constructed pointwise by scanning through 0 ≤ x T ≤ l/2. Results with x P < 0 are discarded. As it can be seen from Figure 9, the algorithm produces a non-zero slope of the contour r (x) in the center of the ligament.
For the stiffness calculation, the ligament is fixed at one end and, due to symmetry, a force F normal to the ligament axis is applied in the symmetry axis, see Appendix of [6] for more details. The bending stiffness is calculated by dividing the applied force F by the deflection w.
Following Euler-Bernoulli beam theory, the second derivative of the deflection w is: which depends on the bending moment M(x) = F(l/2 − x), Young's modulus E and the moment of inertia I(x) = π 4 r(x) 4 . By integrating Equation (9) twice, the deflection w of the ligament is calculated, and thus the bending stiffness k results from Equation (8).
In what follows, a cylindrical reference ligament is chosen having a representative aspect ratio derived from the previous analyses (see Section 3) with the cylinder radius r 0 = 200 nm and ligament length l = 700 nm. The resulting reference volume V 0 = πr 2 0 l is kept constant describing a material with a defined solid fraction but still allowing for different ligament shapes. For the variation of the parabolic ligament geometry, the length remains constant, whereas the a value in Equation (1) was changed stepwise, and the b value was adjusted ensuring a volume equal to V 0 with an error of ±1%, by numerical integration along the whole ligament including the transition into the nodal volume.
Based on the given ligament shape, the appearing shape as determined by the "Thickness" algorithm can be mathematically constructed by Equations (4)- (7), and the geometry parameter r * sym is calculated from Equation (2) by averaging the radius data within the three sections consistently to the approach in Section 3.3. Figure 10 shows the ligament contours for exemplary r * sym values. The results of the volumetric ratio V /V and the stiffness ratio k /k of the apparent shape determined by the "Thickness" algorithm and the original given shape are presented in Table 3. The geometry was changed using values for r * sym in the range from 1.0 (cylindrical) down to 0.7 (concave) including the major part of ligaments, as shown in the histogram in Figure 5a. It should be noted that the "Thickness" algorithm produces overestimated radii also for convex and s-shaped ligaments out of the nature of fitting the biggest sphere. ligament length = 700 nm. The resulting reference volume = is kept constant describing a material with a defined solid fraction but still allowing for different ligament shapes. For the variation of the parabolic ligament geometry, the length remains constant, whereas the value in Equation (1) was changed stepwise, and the value was adjusted ensuring a volume equal to with an error of ±1%, by numerical integration along the whole ligament including the transition into the nodal volume.
Based on the given ligament shape, the appearing shape as determined by the "Thickness" algorithm can be mathematically constructed by Equations (4)- (7), and the geometry parameter * is calculated from Equation (2) by averaging the radius data within the three sections consistently to the approach in Section 3.3. Figure 10 shows the ligament contours for exemplary * values. The results of the volumetric ratio / and the stiffness ratio / of the apparent shape determined by the "Thickness" algorithm and the original given shape are presented in Table 3. The geometry was changed using values for * in the range from 1.0 (cylindrical) down to 0.7 (concave) including the major part of ligaments, as shown in the histogram in Figure 5a. It should be noted that the "Thickness" algorithm produces overestimated radii also for convex and s-shaped ligaments out of the nature of fitting the biggest sphere.
(a) (b) (c) Figure 10. Contours (red curves) obtained from application of the "Thickness" algorithm (red circles), to concave ligaments (blue curve) for geometry parameter values * of (a) 0.9; (b) 0.8; and (c) 0.64. Every 10th calculated construction is marked. All analyzed geometries show surprisingly high stiffness ratios / . Already an * value of 0.91, which is still considered as approximately cylindrical by our 10% tolerance, shows a stiffness overestimation by a factor of 1.43. The / ratio increases gradually with decreasing * values. The maximum stiffness overestimation of a factor of 8.75 occurs at the lowest analyzed * value of 0.64, whereas the average radius of this ligament is overestimated by a factor of 1.32.
In conclusion, these findings underline the importance of the diameter algorithm, especially for the mechanical analysis. Considering just the ligament radii, the "Thickness" algorithm leads to an error of 14% for a concave ligament showing a radius variation of only 10% in * . This has a strong Figure 10. Contours (red curves) obtained from application of the "Thickness" algorithm (red circles), to concave ligaments (blue curve) for geometry parameter values r * sym of (a) 0.9; (b) 0.8; and (c) 0.64. Every 10th calculated construction is marked. Table 3. Results obtained from application of the "Thickness" algorithm to concave ligaments in dependence of the geometry parameter r * sym . Presented are mean radius ratio r / r , volumetric ratio V /V, and the stiffness ratio k /k of the identified to the original ligament shape. All analyzed geometries show surprisingly high stiffness ratios k /k. Already an r * sym value of 0.91, which is still considered as approximately cylindrical by our 10% tolerance, shows a stiffness overestimation by a factor of 1.43. The k /k ratio increases gradually with decreasing r * sym values. The maximum stiffness overestimation of a factor of 8.75 occurs at the lowest analyzed r * sym value of 0.64, whereas the average radius of this ligament is overestimated by a factor of 1.32.
In conclusion, these findings underline the importance of the diameter algorithm, especially for the mechanical analysis. Considering just the ligament radii, the "Thickness" algorithm leads to an error of 14% for a concave ligament showing a radius variation of only 10% in r * sym . This has a strong impact on the mechanical stiffness due to the moment of inertia where the radius enters by a power of four. Similarly, the overestimated radius leads to an overestimation of the mechanical strength.
As the skeleton beam model is built based on the ligament geometry as determined by the "Thickness" algorithm, it is now evident why this model shows the observed degree of overestimation.
It should be noted that the "Thickness" algorithm is a powerful and fast volume-based algorithm. It works well for rod-like structures (e.g., fibers, bones) [42] and is commonly used in image analysis programs, e.g., the open-software program FIJI [35]. The observed overestimation might also not be a problem when studying the self-similarity etc. of structures, or when comparing mean values or distributions [7,15,34]. However, for the mechanical simulation, the correct diameter distribution along the ligament axis is crucial. In the NPG structure, the algorithm reaches its limits due to the strongly varying diameters along the ligaments. A new diameter algorithm is needed for providing the required predictive capability of the skeleton FEM beam model.
Other algorithms are discussed in literature that might be considered alternatively. Approaches like the "nearest neighbor" (closest point-to-point distance between surfaces) or the "surface normal" (measures the diameter as following the normal vector at the surface), or Euclidean distance transformation (EDT), could be tested. Related algorithms are summed up in [42]. In this comparative study, the "Thickness" algorithm also computed the highest diameter values of tibial cartilage thickness. However, the big advantage of this volume-based technique is that no assumptions regarding structure orientation are needed and it is applicable to 3D tomography data while still being fast.
Due to the strongly changing ligament path directions and diameters in NPG, it would be reasonable to use a "skeleton normal" approach, where the medial axis extraction is used for the normal generation [42]. The challenge would be to compute representative normal vectors along the skeleton-line. Such a concept would require a prior smoothening of the skeleton before normal calculation to prevent abrupt changes of the surface search direction along the ligament axis. Here, the treatment of short ligaments is challenging, because only a small amount of voxels are available for generating a smooth skeleton line.
Because of all these complexities, it is preferable to use the "Thickness" algorithm and add a correction approach. As demonstrated in Figure 10, it is technically possible to relate the determined diameter distribution to the underlying original ligament shape. Ongoing work aims at a generalized approach based on inverse methods such as optimization or machine learning that allows tracing back an identified ligament shape to the corresponding true geometry. In view of the diversity of ligament shapes found in Section 3, such an approach should allow for identifying general ligament shapes.

Conclusions
In this paper, a computational efficient skeleton FEM beam model was developed, which is based on skeletonization of 3D FIB-SEM tomography data of NPG. An object-oriented programming approach was chosen allowing for a statistical non-local analysis of the geometrical skeleton network as well as for the generation of a skeleton beam model. The advantage of choosing a modeling approach with beam elements is to save complicated meshing and computational time, along with still having a realistic as possible model by using the tomography data. For comparison, a solid FEM model is created from the very same 3D tomographic data. The computational time is 0.0218 h and 1.38 h for the skeleton beam model and solid model, respectively, both as an elastic-plastic simulation.
A geometrical analysis of the diameter, branch length, and ligament shape was done for the skeleton network produced by the open-source software FIJI. Half of the free-dangling skeleton-branches do not represent significant features of the structure, as they are shorter than the ligament radius. These unwanted branches are produced by the skeletonization algorithm at local ligament thickenings. The remaining longer free-dangling branches are assumed to be pinched-off ligaments formed during the annealing process. Close to the boundary of the RVE in a 3% wide region of the RVE size, a fourfold number of dangling branches is found. Their distances from the boundary and their comparatively small diameter are due to the inclination of ligaments penetrating the RVE plane. Such ligaments need to be included in the boundary conditions although ending slightly below the surface of the RVE.
A high diversity of ligament shapes are found in the NPG structure. The investigation of the ligament shapes revealed a surprisingly high number of s-shape contours: 47% are s-shaped besides 25% concave, 4% convex, and 24% cylindrical ligaments, assuming that deviations below 10% ligament radius are accepted as cylindrical. Because almost half of all analyzed ligaments in the NPG sample are found to have an asymmetric contour, parametrizations and model results limited to symmetrical ligament geometries probably lead to biased statistics in favor of a cylindrical shape.
The comparison of macroscopic stress-strain responses for a compressive strain of 8% revealed that the skeleton beam model overestimates the stiffness by 30% and the plastic strength by 105% in comparison to the solid FEM model. The overestimation of the values by the skeleton beam model was not expected, it is furthermore contradicting the results of the nodal correction proposed by Jiao et al. [17]. This can be traced back to an overestimation of the diameters of the tomography reconstruction by the "Thickness" algorithm for geometries with strongly varying diameters, as they are found in NPG. A study on the influence of the diameter overestimation showed a strong impact on the NPG ligament stiffness. For a typical ligament with a radius variation of up to 0.72 from mid to end radius, a factor of 4.2 in the ligament stiffness was found. It can be concluded that similar factors apply to the macroscopic mechanical properties.
As a next step, a generalized approach based on inverse methods, such as optimization or machine learning, should be developed that allows tracing back an identified ligament shape to the corresponding true geometry. This in combination with the nodal correction is expected to yield much better agreement between the skeleton FEM beam model and the solid FEM model. This will allow producing larger data sets consisting of the complex structure information and related mechanical properties for studying the structure-property relationships of NPG in the future.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
An overview on the analysis (Section 3) and modeling procedure (Section 4) is given in the flow diagram in Figure A1. As input data, various tomography data formats of, e.g., FIB-SEM, synchrotron, or even artificially generated data are possible in general. Figure A1. Flow chart of the analysis and modeling procedure for geometrical and mechanical analysis of NPG networks based on 3D tomography data. For programming the skeleton FEM beam model and the later data analysis, an object-oriented programming (OOP) approach was chosen. This has the advantage of storing the attributes of one object all together with the data processing object functions. Python was chosen as programming language. The scripted Python routine combines the skeleton and diameter data and creates a voxel-object for each voxel of the skeleton with its x-, y-, z-coordinates and the diameter d as attributes.
To model the pairwise relations between the voxel-objects and to achieve the information of the interconnected skeleton network in one skeleton-object, the basics of graph theory are used accordingly for the Python programming, as also in the skeletonization algorithm [43] in FIJI. The skeleton is considered as an undirected graph, which is made up out of vertices interconnected by branches, shown schematically in Figure A2. For establishing the procedure, the tomography slice stack provided by Hu et al. [7,15] was used from which in a first step the skeleton and diameter information are extracted with the open-source program FIJI [35]. The plugin Skeletonize3D [43] included in BoneJ [44] performs the skeletonization based on a thinning algorithm by Lee et al. [37] that converts an image object into its one-voxel-wide inner skeleton centerline. In a second step, the "Thickness" plugin [38] is applied to the original TIFFdata, which constructs the diameter information at each voxel with the "biggest-sphere" algorithm [39].
For programming the skeleton FEM beam model and the later data analysis, an object-oriented programming (OOP) approach was chosen. This has the advantage of storing the attributes of one object all together with the data processing object functions. Python was chosen as programming language. The scripted Python routine combines the skeleton and diameter data and creates a voxelobject for each voxel of the skeleton with its -, -, -coordinates and the diameter as attributes. To model the pairwise relations between the voxel-objects and to achieve the information of the interconnected skeleton network in one skeleton-object, the basics of graph theory are used accordingly for the Python programming, as also in the skeletonization algorithm [43] in FIJI. The skeleton is considered as an undirected graph, which is made up out of vertices interconnected by branches, shown schematically in Figure A2. Figure A2. Schematic 2D drawing of skeleton network and the object-oriented organization of information: The whole skeleton-object is made up out of vertices interconnected by branches. Vertices and branches are made up out of discrete voxels, to which the diameter information is attributed. A branch is always confined by a vertex at both ends, which can be either a junction-or an end-vertex.
A branch is always confined by a vertex on both sides, which can be either a junction, represented in blue, or an end, represented in red. Accordingly, each branch-object knows its two vertex-objects V and V and, beyond classical graph theory, its branch-voxels. A vertex-object knows its junction-voxels or else end-voxel and the projecting branch-objects from this vertex ( , , …, ). To organize a skeleton in the described class structure, each skeleton-voxel is visited and its 3D 26-voxel cubic neighborhood is analyzed, which contains the 26 surrounding voxels with joint faces (6 voxels), edges (12 voxels), and corners (8 voxels). Depending on the number of skeletonvoxels in this neighborhood, the respective type branch-, junction-, or end-voxel is assigned. As in the "AnalyzeSkeleton" algorithm by FIJI [36], end-voxels have less than two neighbors, junctionvoxels have more than two neighbors, and branch-voxels have exactly two neighbors. With this information, the vertex-and branch-objects are created, and thus the initialization of a skeleton-object from the tomography data is completed.
This data structure supports the generation of an ABAQUS model [40] according to the required Figure A2. Schematic 2D drawing of skeleton network and the object-oriented organization of information: The whole skeleton-object is made up out of vertices interconnected by branches. Vertices and branches are made up out of discrete voxels, to which the diameter information is attributed. A branch is always confined by a vertex at both ends, which can be either a junction-or an end-vertex.
A branch is always confined by a vertex on both sides, which can be either a junction, represented in blue, or an end, represented in red. Accordingly, each branch-object knows its two vertex-objects V 1 and V 2 and, beyond classical graph theory, its branch-voxels. A vertex-object knows its junction-voxels or else end-voxel and the projecting branch-objects from this vertex (B 1 , B 2 , . . . , B n ). To organize a skeleton in the described class structure, each skeleton-voxel is visited and its 3D 26-voxel cubic neighborhood is analyzed, which contains the 26 surrounding voxels with joint faces (6 voxels), edges (12 voxels), and corners (8 voxels). Depending on the number of skeleton-voxels in this neighborhood, the respective type branch-, junction-, or end-voxel is assigned. As in the "AnalyzeSkeleton" algorithm by FIJI [36], end-voxels have less than two neighbors, junction-voxels have more than two neighbors, and branch-voxels have exactly two neighbors. With this information, the vertex-and branch-objects are created, and thus the initialization of a skeleton-object from the tomography data is completed.
This data structure supports the generation of an ABAQUS model [40] according to the required input format, which is explained in Section 4. Besides the ABAQUS modeling, several object functions are programmed, that support the analysis of the geometrical skeleton data also non-locally, with regard to the ligament shape along its axis, explained in Section 3.