convolutional, up-convolutional and merge layers are adjusted accordingly. Furthermore, we point out that the cutouts used for training were taken from image data among all seven time steps, but only from the first 200 slices of each image stack; thus, the remaining 220 slices could be used for validation and testing. In order to increase the efficiency of the available training data, we utilized data augmentation (Goodfellow et al., 2016) i.e., during training, pairs of chosen input and corresponding target cutouts were transformed randomly, yet pairwise in the same manner, via rotations/reflections. In this way we increased the number of available input-target pairs, and, additionally, the predictions of the neural network became more stable with respect to rotated images.

As cost function for the training procedure, we chose the binary cross-entropy (negative log-likelihood) function, see Goodfellow et al. (2016). The U-Net's initial kernel weights were drawn from a truncated normal distribution. Then, training of the kernel parameters was performed with the Adam stochastic gradient descent method (Kingma and Ba, 2015), using 50 epochs with 300 steps per epoch and a batch size of 1. These training hyperparameters were manually tuned, while the batch size of 1 was chosen due to memory limitations. The network was implemented using the Keras package in Python, see Chollet (2015) and training was performed on a NVidia GeForce GTX 1080 graphics processing unit (GPU).

After the training procedure, we applied the CNN, denoted by f , to each of the seven available CT images Ct on an Intel Core i5-7600K CPU; that is, we computed predictions for the grain boundary network, Lˆ t , from

$$
\hat{L}\_t = f(\mathbf{C}\_t), \quad \text{for each } t = 0, \dots, 6. \tag{3}
$$

**Figure 8** (middle column) visualizes the outputs Lˆ t of the network in a cross-section (slice 350) that was not used for training. Initial inspection indicates that the predictions of the neural network become less reliable with increasing time or, equivalently, with decreasing visibility of grain boundaries in the CT data. Nevertheless, the predictions are, even for the final time step, reasonably good.

Since the training and application of the 3D U-Net is coupled with high memory usage, we reduced, as already mentioned, the initial resolution of the CT data. Furthermore, **Figure 5** indicates that a 2D U-Net, which can be used at higher resolutions due

FIGURE 13 | (A) 2D cross-section of a CT image containing reconstruction artifacts and (B) the corresponding prediction of the 3D U-Net.

to fewer memory requirements, is capable of detecting grain boundaries from 2D slices—at least for grain boundaries with a good visibility. Therefore, we trained a 2D U-Net using slices, instead of volumetric cutouts. Since the 2D architecture requires less memory than the 3D U-Net, for training we used patches of size 256 × 256 × 1 voxels which were taken from the CT images being downsampled to the resolution of 480 × 480 × 839 voxels instead of 240 × 240 × 420 voxels. In order to allow the 2D U-Net to learn features at a comparable scale as the 3D U-Net, we increased the depth (as defined above) of the 2D U-Net from 4 to 5. After training, the 2D U-Net has been applied slice-by-slice to the seven image stacks, resulting in volumetric grain boundary predictions. Because the 2D U-Net evaluates consecutive slices independently, the network's output can lead to discontinuous grain boundary predictions, see **Figure 9**. To overcome this, we used a 2D U-Net, which was trained with 2D multichannel images with a size of 256 × 256 × 11 voxels. More precisely, it was trained with sets of 11 consecutive CT slices and a ground truth slice corresponding to the 6th input slice. This way, when predicting grain boundaries in consecutive CT slices, the network receives overlapping and correlated information which reduces the discontinuities in the network's output. In order to give the multichannel U-Net additional information for its grain boundary predictions we did not limit it to slice-by-slice predictions in one single axial direction, namely top-to-bottom, of the image stacks. Thus, to obtain the final grain boundary predictions of the multichannel U-Net, the slice-by-slice predictions are computed in three directions (top-to-bottom, left-to-right, front-to-back). For each CT image stack, this results in three grain boundary predictions, which are than averaged resulting in the final volumetric grain boundary predictions.

As of now the procedures described above do not provide a grain-wise segmentation. More precisely, the outputs of the U-Net architectures are 3D images with voxel values in the interval (0, 1). Therefore, the network predictions must be binarized in order to localize the grain boundaries, which, however, do not necessarily enclose grains completely. Therefore, additional

FIGURE 14 | Segmentation results obtained by a 3D U-Net that was trained only with CT/3DXRD data from time step *t* = 6. (A) Kernel density estimation (blue) of relative errors in grain volume. The red curve is the density of relative errors in volume under the condition that the grain is completely visible in the cylindrical sampling window. (B) Kernel density estimation of relative errors in grain volume obtained by the segmentation procedure for each time step *t* = 0, . . . , 6.

image processing steps must be carried out in order to obtain a full segmentation of individual grains.

To that end, we binarize the grain boundary predictions (**Figure 10A**) of the networks Lˆ t with a manually determined global threshold followed by morphological closing (**Figure 10B**). The binarization is followed by a marker-based watershed transformation (**Figure 10C**), which is performed on the (inverted) Euclidean distance transform of the binary images. In order to reduce oversegmentation in the image obtained by the watershed transformation, a final postprocessing step is carried out in which adjacent regions are merged if the overlap between one region and the convex hull of the neighbor is too large (**Figure 10D**). For more details on the marker selection procedure for the watershed transformation or the postprocessing step we refer the reader to Spettl et al. (2015).

In order to quantitatively compare segmentations of the CT images C0, . . . , C6 obtained by the 3D U-Net, 2D U-Net and multichannel U-Net followed by the postprocessing steps described above with segmentations derived from the 3DXRD measurements, we first match grains among these segmentations. More precisely, each grain GXRD ⊂ R 3 observed in a segmentation obtained by 3DXRD microscopy is assigned to a grain Gseg ⊂ R 3 in the corresponding segmentation of the CT image data. We formulated this as a linear assignment problem (Burkard et al., 2012), which minimizes the sum of the volumes of the symmetric differences of matched grains

$$\nu\_{\text{3}}(G\_{\text{XRD}} \Delta \, G\_{\text{seg}}) = \nu\_{\text{3}}(G\_{\text{XRD}} \, \langle \, G\_{\text{seg}} \rangle + \nu\_{\text{3}}(G\_{\text{seg}} \, \langle \, G\_{\text{XRD}} \rangle, \tag{4})$$

where ν3(·) denotes the volume and GXRD 1 Gseg is the symmetric difference given by

$$G\_{\rm XRD} \,\,\Delta \,\, G\_{\rm seg} = \left( G\_{\rm XRD} \,\,\, \left( G\_{\rm seg} \right) \cup \left( G\_{\rm seg} \,\, \left( G\_{\rm XRD} \right) \right) \,. \tag{5}$$

Thus, we will be able to quantitatively compare pairs of matched grains (GXRD,Gseg) which, in turn, allows a comparison of the presented methods.

# 3.3. Results

Even though the CNNs described in section 3.2 do not provide grain-wise segmentation of CT data, they can significantly enhance CT images such that conventional image-processing techniques can be readily used to obtain a grain-wise segmentation. By following the approach described in Spettl et al. (2015), we obtained grain-wise segmentations of the considered data set, despite its rather indistinct grain boundaries. A visual comparison between grain boundaries extracted from the segmentation utilizing a 3D U-Net and the true grain boundaries obtained by 3DXRD microscopy indicates that the segmentation is reasonably good, with some oversegmented grains remaining, see **Figures 10E,F**. A more quantitative comparison becomes available by the grain matching procedure described in section 3.2, i.e., we will compute quantities to measure how much grains segmented from CT deviate from matched grains observed in the ground truth data. More precisely, we determine for pairs of matched grains the relative errors rV in grain volume given by

$$r\_V = \frac{|\nu\_3(G\_{\rm XRD}) - \nu\_3(G\_{\rm seg})|}{\nu\_3(G\_{\rm XRD})}.\tag{6}$$

Also, we computed errors rc in grain barycenter location normalized by the volume-equivalent diameter of the grain GXRD. These values are given by

$$r\_{\mathcal{L}} = \frac{\|\mathcal{L}(\mathcal{G}\_{\text{XRD}}) - \mathcal{L}(\mathcal{G}\_{\text{seg}})\|}{\sqrt[3]{\frac{6}{\pi} \,\nu\_{\text{3}}(\mathcal{G}\_{\text{XRD}})}},\tag{7}$$

where k · k denotes the Euclidean norm and c(GXRD),c(Gseg) are the barycenters of the grains GXRD and Gseg, respectively. **Figure 11** visualizes the quartiles of these relative errors in grain characteristics for the segmentation procedures based on the trained 3D U-Net, 2D U-Net and multichannel U-Net. For reference, we also included results obtained by the conventional segmentation procedure without applying neural networks, which was conceptualized for grain boundaries with good visibility and is described in Spettl et al. (2015). These results indicate that the segmentation procedures based on the U-Net architecture perform better then the conventional method. Among the machine learning approaches, the sliceby-slice approach with the 2D U-Net performs worst with a median value for rV of 0.37. This could be explained by the discontinuities of grain boundary predictions for consecutive slices, see **Figure 9**. By enhancing the slice-by-slice approach with the multichannel U-Net, we achieve a significant drop of this error down to 0.21. The segmentation approach based on the 3D U-Net performs best with a median error of 0.14, because it is able to learn 3D features for characterizing the grain boundary network embedded in the volumetric data.

Kernel density estimations (Botev et al., 2010) of the relative errors forthe 3D U-Net approach are visualized in **Figures 12A,B** (blue curves). Furthermore, **Figures 12C,D** depict these densities for each of the seven observed time steps t = 0, . . . , 6. Note that, as expected, the errors show a tendency to grow with increasing time step. In order to analyze possible edge effects, i.e., a reduced segmentation quality for grains located at the boundary of the cylindrical sampling window, we computed error densities only for grains located in the interior of the sampling window, see **Figures 12A,B**. The plots (red curves) indicate that, indeed, the segmentation procedure based on the 3D U-Net works better for interior grains. This effect can be explained by the information that is missing for grains that are cut off by the boundary of the sampling window.

#### 3.4. Discussion

Although our procedures based on preprocessing with CNNs followed by conventional image processing do not lead to perfect grain segmentations, see **Figure 12**, especially the method utilizing the 3D U-Net delivers relatively good results when considering the nature of the available CT data. Furthermore, the neural network is able to reduce local artifacts, like liquid inclusions in the grain interiors, which cause small areas of high contrast far from grain boundaries, see **Figure 8** (first row). Yet, we warn that the predictions of the trained U-Net are prone to error when there are large-scale image artifacts in the input images, as illustrated in **Figure 13**. One possible way to reduce the effect of such artifacts is to consider a modified architecture of the 3D U-Net, with larger kernels or more pooling layers, such that even larger features can be considered.

Nevertheless, without the machine learning approach, i.e., the preprocessing provided by the 3D U-Net, the segmentation of CT data for later measurement time steps with poorly visible grain boundaries is a complex and time-consuming image processing problem. Still, in the presented procedure, conventional image processing, i.e., binarization and the watershed transform, was necessary to obtain a grainwise segmentation of the considered data. Thus, the segmentation techniques considered in sections 2 and 3 show the flexibility of combining the watershed transform with machine learning techniques either for pre- or postprocessing image data for the purpose of segmenting tomographic image data of functional materials.

Note that, in the 3D U-Net approach, there are some machine learning techniques that could have been adopted to further reduce the need for some of the subsequent image processing steps. For example, the binarization step could be incorporated into the network by using the Heaviside step function as an activation function in the output layer. Morphological operations, like the closing operation utilized in the procedure above, could be implemented by additional convolutional layers with non-trainable kernels followed by thresholding. In this way, the necessary postprocessing steps will be considered during the training procedure of the 3D U-Net. Alternatively, by describing a segmentation with an affinity graph on the voxel grid, it is possible to obtain segmented images as the final output of CNNs, see Turaga et al. (2010). Note that such approaches require cost functions which allow a quantitative comparison between segmentations, see e.g., Briggman et al. (2009) and Liebscher et al. (2015). Furthermore, we point out that there are techniques for obtaining a grain-wise segmentation by fitting mathematical tessellation models to tomographic image data using Bayesian statistics and a Markov chain Monte Carlo approach, see Chiu et al. (2013). In our case, such techniques could be applied directly to tomographic or even to enhanced grain boundary images obtained by the 3D U-Net.

Moreover, we note still another possible application of machine learning methods for the analysis of CT image data. In many applications, "ground truth" measurements are destructive, which means that they can be carried out only for the final time step of a sequence of measurements. This limits the available training data for machine learning techniques.

We simulated such a scenario with our data by using solely the CT image C6 and the 3DXRD data L6 of the last measured time step to train an additional 3D U-Net. Analogously to the procedure described in section 3.2, this network was applied to the entire series of CT measurements. The resulting grain boundary predictions were then segmented using the same image processing steps as described in section 3.2. **Figure 14** indicates that the relative errors of grain volumes are comparable to the errors made when considering every time step during training, see **Figure 12**. This result suggests that a "ground truth" measurement of only the final time step would suffice for training in our scenario. Similarly, machine learning approaches might be interesting for the segmentation and analysis of time-resolved CT data in various applications in which "ground truth" measurements cannot be made during experiments, but only afterwards, in a destructive or time-consuming manner.

# 4. CONCLUSIONS

We gave a short overview of some applications in the field of materials science in which we successfully combined methods of statistical learning, including random forests, feedforward and convolutional neural networks with conventional image processing techniques for segmentation, classification and object detection tasks. More precisely, the methods of sections 2 and 3 utilize machine learning as either a pre- or postprocessing step for the watershed transform to achieve phase-, particle- or grainwise segmentations of tomographic image data from various functional materials—showing how flexible the approach of combining the watershed transform with methods from machine learning is. In particular, we presented such an approach for segmenting CT image data of an Al-5 wt.% Cu alloy with very low volume fraction of liquid between grains. In total, we considered seven CT measurements of the sample, between which were interspersed Ostwald ripening steps. Especially at later times, the aggregation of liquid leads to a decrease in contrast of the image data, i.e., grain boundaries become less distinct in the image data, which makes segmentation by conventional image processing techniques quite difficult and unreliable. Therefore, we employed matching grain boundary images—which had been extracted from the same sample by means of 3DXRD microscopy—as "ground truth" information for training various CNNs: a 2D U-Net which can be applied slice-by-slice to entire image stacks, a multichannel 2D U-Net which considers multiple slices at once for grain boundary prediction in a planar section of the image stack and, finally, a 3D U-Net which was trained with volumetric cutouts at a lower resolution. After the training procedure, the U-Nets were able to enhance the contrast at grain boundaries in the CT data. Especially, the 3D U-Net successfully predicted the locations of many grain boundaries that were either missing from the image data or poorly visible. This shows that machine learning methods can facilitate difficult image processing tasks, provided that "ground truth" data is available, e.g., data obtained via additional measurements or manual image labeling. Since the images output by the convolutional neural networks were not themselves grainwise segmentations, we applied conventional image processing algorithms to the outputs to obtain full segmentations at each considered time step and for each presented method. These were compared quantitatively with "ground truth" segmentations extracted from 3DXRD measurements. The resulting relative errors in grain volume and locations of grain centers of mass indicated that the machine learning-based segmentation

procedures worked reasonably well, particularly for grains that were not cut off by the boundary of the observation window. Finally, we trained an additional 3D U-Net only with CT and 3DXRD data obtained during the final time step. This simulated the common scenario in which a "ground truth" measurement can be performed only at the very end of an experiment. The 3D U-Net trained in this manner was applied as before to the entire CT data set, followed by conventional image processing steps, yielding grain segmentations. Quantitative comparison of the latter to segmentations derived from 3DXRD data indicated that the approach produced good results. Even though a trained neural network does not make 3DXRD measurements obsolete, the procedure presented here can potentially reduce the amount of 3DXRD beam time that is needed for accurate segmentation and microstructural analysis. Likewise, we believe that a similar approach might be particularly beneficial whenever nondestructive CT measurements can be carried out in situ, but "ground truth" information can be acquired only by a destructive measurement technique.

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

OF, MN, LP, and MWe reviewed previous results on machine learning for segmentation of image data. Tomographic image data of the AlCu specimen has been provided by MWa and CK. The network training, segmentation and analysis of AlCu CT image data was performed by OF. All authors discussed the results and contributed to writing of the manuscript. CK and VS designed the research.

#### ACKNOWLEDGMENTS

The financial support of the German Research Foundation (DFG) for funding this research project (SCHM997/23-1) is gratefully acknowledged. The authors thank Murat Cankaya for the processing of image data. In addition, the authors are grateful to the Japan Synchrotron Radiation Research Institute for the allotment of beam time on beamline BL20XU of SPring-8 (Proposal 2015A1580).

lithium-ion batteries under increasing compaction. J. Microsc. 272, 96–110. doi: 10.1111/jmi.12749

solid-liquid mixtures. Acta Mater. 54, 2027–2039. doi: 10.1016/j.actamat.2005. 12.038

**Conflict of Interest Statement:** 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.

Copyright © 2019 Furat, Wang, Neumann, Petrich, Weber, Krill and Schmidt. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Computation of Thickness and Mechanical Properties of Interconnected Structures: Accuracy, Deviations, and Approaches for Correction

#### Claudia Richert 1 \*, Anton Odermatt 1,2 and Norbert Huber 1,2

1 Institute of Materials Research, Materials Mechanics, Helmholtz-Zentrum Geesthacht, Geesthacht, Germany, 2 Institute of Materials Physics and Technology, Hamburg University of Technology, Hamburg, Germany

Identifying local thickness information of fibrous or highly porous structures is challenging. The analysis of tomography data calls for computationally fast, robust, and accurate algorithms. This work systematically investigates systematic errors in the thickness computation and the impact of observed deviations on the predicted mechanical properties using a set of 16 model structures with varying ligament shape and solid fraction. Strongly concave, cylindrical, and convex shaped ligaments organized in a diamond structure are analyzed. The predicted macroscopic mechanical properties represent a highly sensitive measure for systematic errors in the computed geometry. Therefore, the quality of proposed correction methods is assessed via FEM beam models that can be automatically generated from the measured data and allow an efficient prediction of the mechanical properties. The results show that low voxel resolutions can lead to an overprediction of up to 30% in the Young's modulus. A model scanned with a resolution of 200 voxels per unit cell edge (8M voxels) reaches an accuracy of a few percent. Analyzing models of this resolution with the Euclidean distance transformation showed an underprediction of up to 20% for highly concave shapes whereas cylindrical and slightly convex shapes are determined at high accuracy. For the Thickness algorithm, the Young's modulus and yield strength are overpredicted by up to 100% for highly concave ligament shapes. A proposed Smallest Ellipse approach corrects the Thickness data and reduces this error to 20%. It can be used as input for a further robust correction of the Thickness data using an artificial neural network. This approach is highly accurate with remnant errors in the predicted mechanical properties of only a few percent. Furthermore, the data from the FEM beam models are compared to results from FEM solid models providing deeper insights toward further developments on nodal corrections for FEM beam models. As expected, the FEM beam models show an increasing overprediction of the compliance with increasing solid fraction. As an unexpected result, the mechanical strength can however be underpredicted or overpredicted, depending on the ligament shape. Therefore, a nodal correction is needed that solves contradicting tasks in terms of stiffness and strength.

Keywords: tomography, skeletonization, thickness correction, artificial neural network, nanoporous gold, trabecular bone, foams, FEM beam model

Edited by:

Nicola Maria Pugno, University of Trento, Italy

#### Reviewed by:

Ercan Gürses, Middle East Technical University, Turkey Douglas Soares Galvao, Campinas State University, Brazil

> \*Correspondence: Claudia Richert<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Materials Science, a section of the journal Frontiers in Materials

Received: 04 February 2019 Accepted: 28 November 2019 Published: 18 December 2019

#### Citation:

Richert C, Odermatt A and Huber N (2019) Computation of Thickness and Mechanical Properties of Interconnected Structures: Accuracy, Deviations, and Approaches for Correction. Front. Mater. 6:327. doi: 10.3389/fmats.2019.00327

# INTRODUCTION

Lacking a detailed morphological and topological description of the microstructure, the structure-property relationship of open-pore materials, such as metal foams, elastomeric foams, or Nanoporous gold (NPG) is commonly described by the Gibson-Ashby scaling law, in which the solid fraction is the most important parameter characterizing the materials morphology (Gibson and Ashby, 1997; Ashby et al., 2000). During the last two decades, the morphological characterization and prediction of mechanical properties of open-pore materials gained increasing attention, thanks to the improving resolution of X-ray, FIB, and TEM micro-/nanotomography instruments, complemented by advancing image processing algorithms and computational modeling techniques. Tomography and FEM simulations on metal and elastomeric foams date back to Nieh et al. (1998), Nieh et al. (2000), and Kinney et al. (2001). A very detailed analysis of cell volume and strut length distributions, number of faces per cell, junctions coordination number and the shape of the most representative cells was carried out by Dillard et al. (2005) based on a 3D quantitative image analysis of open-cell nickel foams under tension and compression loading using Xray microtomography.

First studies based on NPG were conducted by Rösner et al. (2007) using TEM on dealloyed gold leafs. Hu et al. (2016), Mangipudi et al. (2016), and (Ziehmer et al., 2016) analyzed NPG samples of larger volumes, obtained from focused ion beam (FIB) sectioning and scanning electron microscope (SEM) imaging. By these thorough works, a systematic analysis of the NPG morphology in terms of ligament size distribution and connectivity density has become possible for the first time. Because the ligaments are of nanoscale dimension, these investigations are all based on high-resolution SEM images for which techniques for an automated image processing are an asset. Hu et al. (2016) and Mangipudi et al. (2016) use the 3D Biggest Sphere Thickness algorithm by Hildebrand and Rüegsegger (1997) for the estimating the ligament size distribution of 3D volumes.

For the geometrical description of the ligaments in a NPG network, Pia and Delogu (2015) proposed a parabolic shape with a square cross-section connected in cubic nodes. The parameters for the parabolic shape and their statistical distribution were manually determined from 2D SEM images. Badwe et al. (2017) analyzed 2D SEM images using digital image analysis to obtain ligament size histograms that were fit to the Weibull distribution. To obtain the ligament size distribution, they apply the skeletonization and distance map transformation each onto the original binary SEM image, using the open-source software ImageJ. The multiplication of these two results yields the skeleton ascribed with the according diameter at each skeleton-point. Consistent with the results of Rösner et al. (2007) and Hu et al. (2016), the mean ligament distributions were reported to be nearly self-similar for the examined ligament sizes. Stuckner et al. (2017) present a Python package AQUAMI, which automatically analyzes microstructural features from micrographs. The approach is similar to the approach by Badwe et al. (2017), which was independently published, but has no need for manual calculation in ImageJ. The average diameter and diameter distribution of the morphologies in each phase is calculated using a medial axis transform and a distance transform. McCue et al. (2018) use AQUAMI to data-mine NPG 2D images of 28 published manuscripts, regarding mean ligament diameter, length, and solid phase fraction. They point out the difficulty and resulting systematic discrepancies when comparing results gained by different measuring approaches, ranging from manually measuring the thinnest part of the ligament, to computational estimations. Furthermore, as a minimum criterion for meaningful image analysis, they propose to use images with a minimum resolution of at least 10 pixels per ligament diameter, due to the otherwise reported errors.

In summary, two algorithms are found to be dominantly used in literature to estimate the ligament size distribution: The Thickness algorithm, which is able to analyze 3D volumes and the Euclidean distance transformation (EDT), which is applied for analyzing 2D SEM images by Badwe et al. (2017), Stuckner et al. (2017), and McCue et al. (2018). It calculates at each point of the structure the distance to the nearest background point. The Thickness algorithm by Hildebrand and Rüegsegger (1997) is implemented in image analysis programs, such as the opensoftware program Fiji by Schindelin et al. (2012). It calculates the local thickness at a point as the dimeter of the largest sphere, which is completely inside the structure and which contains the evaluated point. The mean thickness is calculated as the volume weighted average of the local thickness. The algorithm is commonly used to estimate the mean trabecular thickness of trabecular bone (Day et al., 2000; Almhdie-Imjabber et al., 2014), or other bone structures (Witkowska et al., 2014), because it is a powerful and fast volume-based algorithm. In the context of NPG the Thickness algorithm has been applied for analyzing 3D tomography data or voxel models by Hu et al. (2016), Mangipudi et al. (2016), Richert and Huber (2018), and Soyarslan et al. (2018a,b).

By the definition of Hildebrand and Rüegsegger (1997), the biggest sphere at a skeleton point pskel does not need to be centered at this point. Liu et al. (2014) show for an object formed by two overlapping disks of different scales that the Thickness algorithm shows a bias toward the larger disk. They furthermore show that an equivalently working Smallest Sphere approach results in the same artifact, but in the opposite direction. The authors propose the definition of the thickness of a point p as the diameter of the maximum inscribed sphere whose circumference is farthest from p. Furthermore, for the skeleton, the property must be satisfied that the thickness at a skeleton point pskel is the diameter of the biggest sphere centered at pskel. They introduced also a star-line-based algorithm, where the thickness at an axial voxel is defined as the minimum-intercept of a straight line with the boundary. The minimum-intercept length measure is highly robust under small random shifts of axial voxels. One drawback of this thickness computation method lies in the increased computation time needed, because interpolated intensity values at multiple sample points have to be computed on individual starlines for each axial voxel. For more details and other thickness approaches see also the literature cited by Liu et al. (2014). The tendency to overpredict the thickness of structures was also reported by Maier et al. (2017) for cartilage thickness, in comparison to other thickness estimation approaches. Such an overprediction is unproblematic when studying the selfsimilarity of structures, or when comparing mean values or distributions. However, for the prediction of mechanical properties using FEM, the correct diameter distribution along the ligament axis is crucial. Richert and Huber (2018) showed that the Thickness algorithm reaches its limits when being applied to typical shapes of NPG ligaments, due to the strongly varying diameter along the ligament axis. The resulting overestimation in ligament radius up to 30% has a strong impact on the predicted mechanical stiffness, which can deviate by a factor of more than two. In their conclusions, Richert and Huber (2018) mentioned the need for a correction method for tracing back an identified ligament shape to the corresponding true geometry, which could be based on inverse methods, such as optimization or machine learning. This important finding has been ignored by Soyarslan et al. (2018b) who used the diameter information as determined from the Thickness algorithm in their beam-FE model, without any local validation of the detected diameters or discussion of possible consequences for their mechanical prediction.

Further literature research revealed that there exists also a plugin in the open-software program Fiji of the 3D Euclidean distance transformation (EDT) by Ollion et al. (2013), among others, which seems to be unnoticed by groups working on the analysis of 3D data. As this algorithm computes the distance from a given voxel of the structure to the nearest background voxel, the extracted axis-to-surface distance will have the tendency to underpredict the ligament diameter for highly convex or concave ligament shapes. The reason for this is that the smallest distance is determined by the normal from the surface contour to an axis point, which is smaller compared to the diameter measured normal to the ligament axis. It his however unclear, how large the deviations are for the typical geometries found in open pore materials and how big their impact is on the mechanical properties in comparison to the results from the Thickness algorithm.

Motivated by these findings, this paper aims to lay a solid basis for error estimation and thickness correction for the different algorithms. The availability of a method for an accurate characterization represents a key element for producing data sets of high quality, consisting of pairs of structure information and related mechanical properties. As demonstrated by Huber (2018) for the topology term of the structure-property relationship, a larger number of such patterns is needed for deriving a fairly general representation using data mining and machine learning approaches. This is particularly an issue when pooling data from different sources, which make use of different algorithms.

Following a detailed investigation of the sources of over- and underestimation in the computed thickness data, approaches for the correction of data from the Thickness algorithm are proposed: A Smallest Ellipse algorithm, which resides in between the Biggest Sphere approach and the Smallest Sphere approach, and an artificial neural network approach. Similarly, an artificial neural network approach is proposed for the correction of data from the Euclidean distance transformation. The results clearly show that the artificial neural network is able to correct the overand underpredicted thickness dependent on the position of the ligament axis. The drawback is that it is limited to the range of ligament shapes used during training. Recommendations are given in terms of generalization to asymmetric ligaments as a requirement for applications to larger structures of higher complexity.

# METHODOLOGY

Previous analysis by Richert and Huber (2018) on actual NPG tomography data produced by Hu et al. (2016) revealed a diameter overestimation of the NPG structure by the Biggest Sphere Thickness algorithm by Hildebrand and Rüegsegger (1997), implemented in the open-source program Fiji by Schindelin et al. (2012), in the Thickness Plugin by Dougherty and Kunzelmann (2007). Richert and Huber (2018) mathematically calculated the influence on the overestimated ligament diameters on the mechanical stiffness for single parabolic ligaments, showing an overestimation by up to a factor of 8. These results clearly show the significance of the error to be expected as function of the ligament geometry, but it is unclear how strong this effect is reflected in the macroscopic properties of a Representative Volume Element (RVE). It can be argued that the macroscopic response of an interconnected structure could be less sensitive to local deviations in the ligament geometries. Furthermore, the amount and effect of possible underestimations by the distance transformation need to be investigated. An impression of the discrepancy between the two algorithms is obtained by analyzing the tomography data of Hu et al. (2016), shown in **Figure 1**. The Thickness (Th) and Euclidean distance transformation (EDT) information are consistently evaluated along the skeleton voxels. It can be seen that the determined averages of 400 nm (Th) and 308 nm (EDT) deviate significantly. It is therefore important to investigate each algorithm with

FIGURE 1 | Ligament diameter distribution of NPG tomography with Thickness (Th) and Euclidean distance transformation (EDT) algorithm. The histograms are normalized to an area of one and fitted with the Gaussian distribution. Shifted distributions with average ligament diameter of 400 nm (Th) and 308 nm (EDT) are observed.

respect to ligament shape and to propose correction methods, where needed.

It should be noted that working with tomography data, several crucial image-processing steps are necessary beforehand, such as image noise filtering, brightness and contrast adjustment, registration and segmentation. For the latter, it is necessary to set a threshold value that decides if a voxel is attributed to the solid or to the pore space and the proper choice of this parameter is absolutely critical for all following steps. Commonly, this parameter is calibrated via the relative density of the material, which is independently measured. While this ensures that the tomography reflects the relative density of the material in average, this does not guarantee that local features are precisely detected. In case of the NPG-epoxy composite tomography data produced by Hu et al. (2016), specific settings in the FIB-SEM process made the ligaments easily distinguishable without interfering with the ligament network structure underneath the cross-section. In this case, the segmentation in Fiji using a single value grayscale threshold for the image stack was thus applicable. An image processing error of ±2% in volume fraction was found by manually changing the image contrast, brightness and threshold value for the segmentation process for that data set (Hu, 2017).

This study focuses on analyzing the influence of the Thickness and EDT algorithm on NPG-like RVEs, which are based on known geometries. Emphasis is placed on providing data of sufficiently complex but well-defined 3D structures, for which the exact diameter information is known in each position along the ligament axis. To this end, ligaments with a smooth parabolic-spherical ligament shape as suggested by Richert and Huber (2018) are organized in a diamond structure. This topology is frequently used for mechanical modeling of 3D open pore materials (Nachtrab et al., 2011; Huber et al., 2014; Roschning and Huber, 2016; Jiao and Huber, 2017a,b; Huber, 2018). In contrast to the conventional FEM approaches, which are computational expensive, FEM beam models allow for fast computation even for large plastic deformation, which is a requirement for larger parameter studies of larger and more realistic RVEs. The drawback of this method is the underprediction of stiffness and strength, which needs to be compensated via a correction of the nodal mass (Huber et al., 2014; Roschning and Huber, 2016; Jiao and Huber, 2017b). An attractive alternative for the numerical simulation of foam-like materials is the Finite Cell Method (Parvizian et al., 2007; Düster et al., 2008, 2017). Recently, Gnegel et al. (2019) applied this approach for predicting the elastic-plastic deformation behavior of pure and polymer coated NPG based on the tomography data of Hu et al. (2016). In combination with experimental macroscopic compression data, it was possible to determine the elastic-plastic properties of the gold phase and of the polypyrrole coating of a few nanometer thickness. This requires reducing the explicitly modeled 3D structure to a sub-sample of the available tomography dataset such that the model could be computed in a reasonable time. Therefore, FEM beam models remain an attractive candidate for computing larger models.

For the sake of a systematic in-depth comparison of all methods under investigation, the geometries in this work are limited to symmetric shapes. Altogether, 16 idealized model geometries plus three additional validation geometries are generated covering the relevant range of ligament shapes from concave to convex. For each model geometry, a high-resolution voxel representation serves as basis for testing various approaches of thickness detection and correction. In addition to the assessment of the error in the determined geometry, the effect on the mechanical properties is computed for each structure and correction method using the FEM beam modeling approach developed in a series of previous works (Huber et al., 2014; Jiao and Huber, 2017a; Huber, 2018; Richert and Huber, 2018).

Motivated by the reported differences between the skeleton FEM beam model and the FEM solid model (Richert and Huber, 2018), FEM solid models are created via PCL scripting in MSC Patran, complementing the reference FEM beam models. The results will provide further insights into the differences between FEM beam and FEM solid models for various ligament shapes in terms of elastic and plastic deformation behavior. The results are also relevant for the further development of nodal corrections for more general ligament shapes as an extension to the simple ball-and-stick geometries investigated by Jiao and Huber (2017b).

**Figure 2** gives an overview of the workflow applied in the following sections. Details on the individual approaches are provided at the beginning of each section. To mimic the FEM skeleton beam model building process from tomography data by Richert and Huber (2018), the RVE geometry information is scanned by a Python script with a defined voxel resolution. The output is a voxelized tiff stack, which is needed as input for the Skeletonize, AnalyzeSkeleton, Thickness and 3D Distance Map Plugin evaluations in Fiji (Lee et al., 1994; Dougherty and Kunzelmann, 2007; Arganda-Carreras et al., 2010; Ollion et al., 2013). The whole procedure of building the FEM skeleton beam model from tomography data is described in detail in the Appendix of Richert and Huber (2018). The simulation of the original FEM beam model vs. the FEM skeleton beam model created with the Thickness information will reveal the impact of the flawed diameter estimation on the mechanical behavior of the ligament network. This allows us also to individually analyze the errors originating from the voxel resolution, the skeletonization, and the ligament discretization on the macroscopic elastic-plastic response.

After the analysis of the influencing parameters with regard to their effect on the geometry computation, the question arises, to what extend the error of each algorithm could be reduced in the aftermath. Concerning the Thickness algorithm we focus in this work on two different correction approaches. Geometrically it is clear why the Thickness algorithm overestimates the diameters of strongly varying ligament shapes as found in NPG. This is why a direct reconstruction approach is developed, opting for an ellipse as the final scanning volume. This so-called Smallest Ellipse (SE) algorithm resides in between the Biggest Sphere approach and the Smallest Sphere approach and is therefore a promising technique for efficiently balancing the thickness data between over and underprediction. A second correction approach is based on an artificial neural network (ANN), which efficiently allows for a global mapping from the measured overpredicted to the corrected ligament shapes. The ANN approach is also applied for correcting data from the EDT algorithm.

and Euclidean distance transformation are done in Fiji. 3rd step: FEM skeleton beam models are built via python scripting with Thickness (Th), Euclidean distance transformation (EDT), and corrected diameters using the Smallest Ellipse (SE). 4th step: additional artificial neural network (ANN) correction approach.

#### REFERENCE FEM MODELS AND THEIR PROPERTIES

#### Reference Geometry of the Unit Cell

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 sshape 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 G11 and G14 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 rend,l and rend,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 xQ,l and xQ,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).

$$r(\mathbf{x}) = \begin{cases} \sqrt{r\_{end,l}^2 - (l/2 + \mathbf{x})^2} & -l/2 \le \mathbf{x} < \mathbf{x}\_{Q,l} \\\ a\mathbf{x}^2 + b\mathbf{x} + \mathbf{c} & \mathbf{x}\_{Q,l} \le \mathbf{x} \le \mathbf{x}\_{Q,r} \\\ \sqrt{r\_{end,r}^2 - (l/2 - \mathbf{x})^2} & \mathbf{x}\_{Q,r} < \mathbf{x} \le l/2 \end{cases} \tag{1}$$

The axial coordinate x has its origin in the mid of the ligament, such that the ligament mid radius is given by rmid = c. For the indepth 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,rend = rend,l = rend,r and xQ,l = − xQ,r .

In what follows, the unit cell size aUC 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 rmid/rend 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 rend = 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 rmid and rend, the parameter c in Equation (1) can be determined following Richert and Huber (2018).

#### Reference FEM Solid and Beam Models

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 (G11) and highest solid fraction (G44), 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 rend = 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 Es = 80 GPa, a Poisson's ratio of ν = 0.42, a yield strength of σy,s = 500 MPa, and a work-hardening rate of ET = 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.

#### Boundary Conditions

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

Two digits numbering the row and column in this table are used for coding the geometry.

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 elasticplastic 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 G11 and G14 (rend = 0.1, rmid/rend = 0.5 and rmid/rend = 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 G14, as can be seen from **Figures 4C,D**. However, for structure G11 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.

# Reference Macroscopic Mechanical Properties

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

structure G11; (C) solid model of structure G14; (D) beam model of structure G14.

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) , σ (ref) y ) 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 rmid/rend, 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% (rmid/rend = 0.5) or 80% (rmid/rend = 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.

# FEM SKELETON BEAM MODELS

The FEM skeleton beam model building approach of Richert and Huber (2018) is based on tomography data sets of real NPG provided by Hu et al. (2016). The common problem for this and similar works (Mangipudi et al., 2016; Soyarslan et al., 2018b) is that the desired thickness information normal to the ligament axis is not easily available. The 16 model geometries, defined in section Reference Geometry of the Unit Cell, enable us to systematically study the different sources of over- and underprediction and to qualify proposed correction methods. Furthermore, the sensitivity with respect to the voxel resolution, the skeletonization, and the discretization of the ligaments is studied.

#### RVE Size and Voxelization

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 scanbox edge length around the mid-point is limited to 1.5 times of the unit cell size aUC, 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 G11 with the smallest diameter was chosen to determine the accuracy as function of the voxel resolution. This most filigree structure with rend = 0.1 and rmid/rend = 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 Nv/aUC (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 aUC = 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**).

# Skeletonization and Beam Discretization

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

ligament (Bez 20 E/L). The values are fit to a simple hyperbolic function E (Th)= k (Nv/aUC) +E (Th) ∞ . The parameter E (Th) ∞ approximates the values for an infinite number of voxels Nv/aUC, being 590 and 580 MPa for Bezier and 1 V/E discretization, respectively. The percentage deviation from those values is inscribed.

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 G11 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 Nv/aUC = ∞, 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 Nv/aUC = ∞. 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 Nv/aUC = ∞. 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 Nv/aUC 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.

#### Thickness Analysis

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 D r (.)E obtained from averaging over all 20 elements of a ligament normalized by the mean-radius of the reference geometry D r (ref) E . It can be seen that the deviation of D r (.)E / D r (ref) E increases with increasing concavity, independent of the end radius rend 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 G41, the deviation is −30%.

For the Thickness algorithm, the local overestimation of the rmid 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 Rnode ≥ rend, forming a nodal plateau of constant radius. Secondly, from there the ligament shape assumes a smooth transition from Rnode to rmid. 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 r (Th) mid can take all possible values from r (ref) mid to Rnode.

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.

#### Effect on Mechanical Properties

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 rmid/rend, independent of the absolute rend 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 rend values are observed, especially for the concave ligaments. There, smaller rend 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

FIGURE 10 | Results of the macroscopic mechanical properties for the FEM skeleton beam models based on the Thickness (Th) algorithm or Euclidean distance transform (EDT), normalized to the results from the reference FEM beam models: (A) Young's modulus E and (B) yield strength σy . The superscript (.) corresponds to (Th) or (EDT).

moment, which again depends on the lever acting on this crosssection. 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.

# METHODS FOR THICKNESS CORRECTION

Due to the impact on the mechanical response, we present in the following sections possible correction approaches for both thickness algorithms. The high sensitivity of the mechanical properties on the geometric characterization justifies to use the predicted Young's modulus and yield strength throughout these sections as the relevant measure for the assessment of the quality of each approach.

# Smallest Ellipse Approach

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.

# Artificial Neural Network Correction Approach

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

Young's modulus; (B) Macroscopic yield strength at 1% plastic strain.

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)/ D r (.)E , 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 2xi/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 xi . 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.

$$\mathbf{X} = \left\{ \frac{r^{(.)}\left(\mathbf{x}\_{1}\right)}{\left\langle r^{(.)}\right\rangle}, \frac{r^{(.)}\left(\mathbf{x}\_{2}\right)}{\left\langle r^{(.)}\right\rangle}, \dots, \frac{r^{(.)}\left(\mathbf{x}\_{20}\right)}{\left\langle r^{(.)}\right\rangle}, \frac{2\mathbf{x}\_{i}}{l} \right\} \tag{2}$$

$$\mathbf{Y} = \left\{ \frac{r^{(ref)} \left( \mathbf{x}\_i \right)}{r^{(\cdot)} \left( \mathbf{x}\_i \right)} \right\} \tag{3}$$

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 rend = 0.1375 and rmid/rend = [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 r ∗ sym and r ∗ asym, 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 underpredictions 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.

# APPLICATION TO EXPERIMENTAL TOMOGRAPHY DATA

To test the methods presented in section Methods for Thickness Correction beyond the 16 idealized diamond structures, the NPG tomography data set of Hu et al. (2016) is used, which stems from a nanoporous gold sample with nominally 400 nm average ligament diameter. Three FEM skeleton beam models are generated based on the Thickness, EDT and furthermore Smallest Ellipse corrected data of the NPG tomography, as described in section FEM Skeleton Beam Models. The ANN approach is not applicable in its present form, as it requires an extension toward more general shapes. The mesh of the reference Solid model consisting of 10-node tetrahedral elements (C3D10) was provided by Hu et al. (2016). For both types of models, symmetry boundary conditions are used and a compressive loading in zdirection is applied. In this way, the results give an additional insight about the effect for a more commonly used boundary condition and a realistic, aperiodic microstructure. The resulting macroscopic stress-strain behaviors are plotted in **Figure 14**; the inserts clearly show the differences in the resulting beam models relative to the solid model, where the black line shows the traced outline of the solid model.

The macroscopic Young's moduli and yield strengths at 1% plastic strain are computed as 432 and 10.0 MPa (Thickness); 310 and 7.5 MPa (Smallest Ellipse); 160 and 3.7 MPa (EDT), respectively. For comparison, the Young's modulus and yield strength of the FEM solid model was computed as 370 and 5.5 MPa. Consistent with the trend observed for the idealized diamond structures, the model based on Thickness and EDT diameter information show the highest and lowest values, respectively. From the idealized diamond structures, we computed ratios (Thickness/EDT) in the average diameter ranging from 1.05 (highest convexity) to 1.4 (highest concavity). The corresponding ratios in the predicted mechanical properties are ranging from 1.2 to 2.5 (Young's modulus) and 1.2 to 2.7 (yield strength). For the NPG tomography, the ligament diameter distribution resulting from the Thickness and EDT algorithms showed a ratio in the average ligament diameter of 1.3 (see section Methodology and **Figure 1**). As shown in **Figure 14**, the (Thickness/EDT) ratio computed for the macroscopic mechanical properties are 2.7 (Young's modulus) and 2.7 (yield strength). Therefore, geometry and property ratios for the real material are close to or above the upper limits found for the idealized diamond structures. This is reasonable, because the diamond structure exhibits straight skeleton lines and symmetric ligament profiles while the skeleton paths in the NPG sample are more randomized and ligament profiles are strongly asymmetric, as reported by Richert and Huber (2018). Therefore, the ligaments show additional gradients along their axis—such gradients are found to be the source of error in both algorithms, Thickness and EDT.

Furthermore, the resulting Young's modulus of the EDT beam model is only 43% of the solid model. This confirms the expectation by Richert and Huber (2018) that the FEM Solid model should be stiffer and stronger than the FEM skeleton beam model, as in the latter, the stiffening and strengthening effect of the nodal mass (Jiao and Huber, 2017b) is not yet accounted for.

# CONCLUSIONS AND OUTLOOK

While the accurate determination of the thickness of geometrical features from 2D images is straight forward, the situation changes dramatically for 3D structures. Various algorithms exist, but each has its specific drawbacks regarding implementation, computational cost, or accuracy. The Thickness algorithm by Hildebrand and Rüegsegger (1997) is the most commonly used algorithm. This is usually done without an assessment of the error, because information about the correct thickness of the structures under investigation is not available. A study by Richert and Huber (2018) of typical ligament shapes identified from 3D FIB tomography data of NPG revealed that the error in the measured geometry can reach values up to 30%. The overestimated thickness data lead to an overestimation of the mechanical stiffness by a factor of two and more. Although an implementation of the 3D Euclidean distance transformation (EDT) is for example available in the Plugin TANGO, this algorithm has so far not been used in 3D analysis. In contrast to the Thickness algorithm, it tends to underpredict the diameter for curved shapes. A first comparison of both algorithms with tomography data of NPG revealed a difference in the computed average ligament diameter of 30%. This and the detailed results obtained on the local radii for the different algorithms highlight how important it is to understand the individual algorithm used and what the produced data represent in relation to the measure of interest. This is particularly an issue when pooling data from various sources making use of different algorithms.

To provide RVEs with well-defined geometries, this work is based on idealized model structures consisting of ligaments with circular cross-sections and smooth parabolic-spherical shape, organized in a diamond structure. Sixteen high-resolution voxel models and finite element models are provided, covering the relevant shapes from concave to convex ligaments. These models serve as reference for the error assessment for both, the determined geometry and the elastic-plastic mechanical properties along the thickness determination and correction chain. Furthermore, the provided test structures can be used for validation of any newly developed algorithm for the determination or correction of thickness information from voxel data.

FIGURE 14 | Stress strain curves of NPG tomography of Hu et al. (2016) modeled as FEM solid model and predictions from FEM skeleton beam models based on Thickness, Smallest Ellipse, and EDT diameter information. Inserts show regions of interest for comparison of the determined geometries.

To decouple this study from the known effect that FEM beam models show a more compliant and less strong macroscopic stress-strain behavior compared to FEM solid models, the differences in both properties are computed for each geometry. As expected, the FEM beam model is more compliant compared to the FEM solid model. The data show an increasing deviation for increasing mid to end radius ratio while the ligament size has only a marginal effect. In contrast to this, the yield strength distributes below and above those of the FEM solid model. This surprising result leads to the conclusion that the stressstrain curve computed by Richert and Huber (2018) must not necessarily fall below the curve predicted by the FEM solid model, after the geometry is corrected, because a newly developed nodal correction for these ligament shapes may not necessarily increase the strength.

An investigation of the sensitivity with regard to the voxel resolution revealed that the predicted mechanical stiffness is significantly overestimated with decreasing voxel number. For the most filigree structures and a resolution of 60 voxels per unit cell length, the error reaches up to 30% in comparison to a resolution of 300 voxels. Increasing the resolution to 200 voxels reduces the error to 3%.

Applying the Thickness algorithm to the data with 200 voxels resolution yielded largest overestimations of 20% in the average radius and 70% in the radius in the middle of the ligament. The impact on the Young's modulus and yield strength is up to 100% overestimation for concave shapes. This is not as high as predicted in the single ligament study by Richert and Huber (2018), but is still inacceptable. The Euclidean distance transformation resulted in an underprediction in the macroscopic mechanical properties of up to 20% for concave ligaments.

In view of these results, two approaches for correction of the computed thickness are proposed. A Smallest Ellipse correction approach, which could be interpreted as counterpart of the Biggest Sphere Thickness algorithm, allows reducing the error in Young's modulus to 20% and in yield strength to 30% for all ligament shapes. Secondly, using patterns consisting of estimated thickness information from Thickness, Smallest Ellipse, Euclidean distance transformation algorithm, and original ligament shapes, artificial neural networks were trained. It could be shown that the accuracies achieved for most cases are within a few voxels. The resulting deviations in the mechanical properties are within few percent, even for untrained validation patterns. This demonstrates the big potential of ANNs to accurately approximate complex non-linear relationships as whole. Even a correct reconstruction is possible for data for which the input information is incomplete in terms of the original ligament shape. However, relative to the ANN corrected Thickness data, the accuracy can be significantly increased by presenting the data from the Smallest Ellipse algorithm. This shows that it is advisable to reduce the complexity of the problem as far as possible by using existing algorithms or estimates, even if they are of limited accuracy. Such strategies have been successfully applied before and the outcome of this work emphasizes once more the importance of incorporation of a priori knowledge in the preparation of the ANN definition and pattern generation when high accuracy is a requirement. This is particularly important for solving highly non-linear and complex inverse problems (Huber and Tsakmakis, 1999, 2001; Tyulyukovskiy and Huber, 2006).

An obvious drawback of the ANN approach is that it must be trained for the parameter space of possible shapes to be identified. This means that for the evaluation of tomography data in the parameter space of r ∗ sym and r ∗ asym, as found by Richert and Huber (2018), requires an expansion by a linear gradient along the ligament axis or the incorporation of even more general shapes. In addition, the results of the Thickness and EDT algorithms should be critically evaluated with respect to effects from non-circular cross-sections that might occur in real samples.

Thus, future research should be directed toward approaches that provide sufficient geometrical accuracy for a large range of possible ligament geometries, where the accuracy should always be evaluated in view of the predicted mechanical properties.

Finally, the Thickness, Smallest Ellipse, and EDT algorithms are applied to the experimental NPG tomography data set of Hu et al. (2016). The average diameters and predicted stressstrain curves consistently showed Thickness to EDT ratios at the upper limit of the range computed for the idealized diamond structures. This is consistent with the finding that gradients in the ligament diameter along the axis are responsible for systematic over- and underestimation by the algorithms. Obviously, this effect is enhanced by the random nature and strong asymmetry of real ligaments. Furthermore, the stressstrain curve of the solid model lies in between the Thickness and EDT prediction. While the overprediction based on Thickness data confirms the result reported by Richert and Huber (2018), the EDT curve being significantly below the result of the FEM solid model now opens the perspective for an implementation of a physically meaningful nodal correction in the FEM beam model.

# AUTHOR CONTRIBUTIONS

CR and NH conceptualized, designed the study, and wrote and revised the manuscript. NH created the geometries, FEM solid and FEM beam models, and coded the Python scripts for voxel scanning. CR analyzed the voxel scans, carried out the FEM skeleton beam models. AO developed the Smallest Ellipse approach and carried out the thickness correction using this method. CR developed the ANN for thickness correction, analyzed the errors of all methods in terms of geometry and mechanical properties.

# FUNDING

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Projektnummer 192346071—SFB 986.

# ACKNOWLEDGMENTS

Kaixiong Hu and Erica T. Lilleodden are acknowledged for making the FIB tomography data set and the FEM solid model of NPG available.

### REFERENCES

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmats. 2019.00327/full#supplementary-material

quantitative analysis of morphologically complex multiphase materials. Comput. Mater. Sci. 139, 320–329. doi: 10.1016/j.commatsci.2017. 08.012

**Conflict of Interest:** 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.

Copyright © 2019 Richert, Odermatt and Huber. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Data Science Based Mg Corrosion Engineering

Tim Würger 1,2, Christian Feiler 1 , Félix Musil 3 , Gregor B. V. Feldbauer 4 , Daniel Höche1,5 , Sviatlana V. Lamaka1 , Mikhail L. Zheludkevich1,6 and Robert H. Meißner 1,2 \*

1 MagIC - Magnesium Innovation Centre, Institute of Materials Research, Helmholtz Centre for Materials and Coastal Research, Geesthacht, Germany, 2 Institute of Polymers and Composites, Hamburg University of Technology, Hamburg, Germany, 3 Laboratory of Computational Science and Modelling, Institute of Materials, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, 4 Institute of Advanced Ceramics, Hamburg University of Technology, Hamburg, Germany, 5 Computational Material Design, Faculty of Mechanical Engineering, Helmut-Schmidt-University, Hamburg, Germany, 6 Faculty of Engineering, Institute for Materials Science, University of Kiel, Kiel, Germany

Magnesium exhibits a high potential for a variety of applications in areas such as

#### Edited by:

Surya R. Kalidindi, Georgia Institute of Technology, United States

#### Reviewed by:

Tim Kowalczyk, Western Washington University, United States Ruifeng Lu, Nanjing University of Science and Technology, China

> \*Correspondence: Robert H. Meißner<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Materials Science, a section of the journal Frontiers in Materials

> Received: 04 January 2019 Accepted: 14 March 2019 Published: 05 April 2019

#### Citation:

Würger T, Feiler C, Musil F, Feldbauer GBV, Höche D, Lamaka SV, Zheludkevich ML and Meißner RH (2019) Data Science Based Mg Corrosion Engineering. Front. Mater. 6:53. doi: 10.3389/fmats.2019.00053 transport, energy and medicine. However, untreated magnesium alloys are prone to corrosion, restricting their practical application. Therefore, it is necessary to develop new approaches that can prevent or control corrosion and degradation processes in order to adapt to the specific needs of the application. One potential solution is using corrosion inhibitors which are capable of drastically reducing the degradation rate as a result of interactions with the metal surface or components of the corrosive medium. As the sheer number of potential dissolution modulators makes it impossible to obtain a detailed atomistic understanding of the inhibition mechanisms for each additive, other measures for inhibition prediction are required. For this purpose, a concept is presented that combines corrosion experiments, machine learning, data mining, density functional theory calculations and molecular dynamics to estimate corrosion inhibition properties of still untested molecules. Concomitantly, this approach will provide a deeper understanding of the fundamental mechanisms behind the prevention of corrosion events in magnesium-based materials and enables more accurate continuum corrosion simulations. The presented concept facilitates the search for molecules with a positive or negative effect on the inhibition efficiency and could thus significantly contribute to the better control of magnesium / electrolyte interface properties.

Keywords: machine learning, property-structure relationship, high-throughput screening, corrosion inhibition, density functional theory, magnesium, dimensionality reduction

# 1. INTRODUCTION

Light-weight materials such as magnesium and its alloys are of high interest for the industrial sector. Potential applications can be found in the automobile industry as structural component (Kulekci, 2008), in batteries as anode material (Aurbach et al., 2000; Höche et al., 2018) and in medical engineering as biocompatible, resolvable implant (Brar et al., 2009). However, dealing with corrosion is a challenging task in various engineering disciplines. Durability and versatility strongly depend on the corrosion properties of the applied material and for most applications as structural component, corrosion activity has to be minimized. Yet, for other approaches the corrosion properties have to be adapted to fit the desired application. As for example, introducing magnesium as battery or implant material requires the corrosion or degradation to proceed with a certain rate. Consequently, the development of reliable, predictive models and methods for general dissolution control is crucial.

There are several concepts to protect magnesium from corrosion, ranging from alloying to surface coatings (Gray and Luan, 2002; Blawert et al., 2006; Jia et al., 2016). Recent studies strongly suggest that the re-deposition of released noble impurities (e.g., iron) results in higher corrosion rates as the size of cathodically active sites at the magnesium surface increases over time up to a state of equilibrium (Höche et al., 2016; Li et al., 2016; Mercier et al., 2018; Michailidou et al., 2018). Concerning the iron re-deposition mechanism, a promising strategy to prevent or control corrosion in magnesium-based materials is the introduction of chemical substances that either form stable complexes with the released iron species or block their access to the surface (Lamaka et al., 2016; Yang et al., 2018).

Novel methods for inhibition prediction of not yet tested compounds based on modern data science techniques are in high demand to predict whether a molecule is a potential inhibitor or even further promotes dissolution of the used material. Hence, the high experimental effort and costs of testing multiple compounds for their corrosion inhibition potential can be circumvented. The molecular structure of potential corrosion inhibiting additives is easily obtained nowadays and thus, represents a promising starting point to identify property-structure relationships as well as to predict the inhibition efficiencies of uninvestigated additives. Following this strategy, Ceriotti et al. developed sophisticated methods to vividly illustrate property-structure landscapes by employing SOAP (Smooth Overlap of Atomic Positions) kernels (Bartók et al., 2013) to create a high-dimensional similarity measure and reducing it to a two-dimensional visualization with the dimensionality reduction algorithm "sketch-map" (Ceriotti et al., 2011, 2013). Moreover, this approach is particularly suited for high-dimensionality data from atomistic simulations as it was already successfully applied to molecular crystals (Musil et al., 2018) and high-throughput structural databases (De et al., 2016, 2017).

In this study, the capabilities of the SOAP kernel and sketch-map are focused on a corrosion inhibition database for multiple molecular compounds to improve the understanding of the inhibition-structure relationship. Furthermore, obtained results can be directly used to qualitatively predict the inhibition properties of not yet tested compounds, thus allowing for a data-driven design of anti-corrosion additives for magnesium-based materials.

# 2. MATERIALS AND METHODS

#### 2.1. Corrosion Experiments

The balance between magnesium dissolution and hydrogen evolution dominates the aqueous magnesium corrosion process. Due to the processing of magnesium with various methods (Pekguleryuz et al., 2013), noble impurities, as for example iron, are impossible to avoid. Thus, local galvanic cells are induced into the material that locally promote the corrosion, resulting in increased magnesium dissolution, hydrogen evolution and the release of impurities, such as iron. Finding molecules that form stable soluble or insoluble complexes with the released impurities is a promising way to screen for dissolution modulators and provides the basis for our workflow.

In a systematic screening for magnesium corrosion inhibitors (Lamaka et al., 2017), the influence of various organic molecules on the hydrogen evolution rate in magnesium corrosion was investigated. Here, the compounds were either previously reported as magnesium corrosion inhibitors or chosen based on their ability to form stable soluble complexes with Fe2+/3+ in order to prevent iron re-deposition (Höche et al., 2016; Lamaka et al., 2016). Hydrogen evolution tests were performed for six different alloys as well as three grades of pure magnesium. Based on the resulting hydrogen evolution rate, the inhibitors were ranked by their inhibition efficiencies, where positive values, up to 99% correspond to suppressed Mg corrosion (referred to as corrosion inhibitors) and negative values to promoted dissolution of Mg (referred to as corrosion promoters) with respect to a reference experiment in 0.5% NaCl electrolyte without any additives. The potential inhibitors were dissolved in 0.5% NaCl to obtain concentrations of 0.05 M and the initial pH was adjusted to the values in the range of 5.5 − 7.2. Further experimental details can be found in the original publication (Lamaka et al., 2017). In this study, only inhibition results for commercial purity magnesium (CP-Mg) with 220 ppm iron content are considered.

# 2.2. Molecular Similarity

SMILES (simplified molecular-input line-entry system) strings of the experimentally investigated compounds are used to create molecular structures using the small molecule topology generator STaGE (Lundborg and Lindahl, 2015). As implemented in the high-throughput workflow of STaGE, the structures are geometry optimized with GAMESS/US (Schmidt et al., 1993; Gordon and Schmidt, 2005) using the B3LYP functional (Becke, 1993; Stephens et al., 1994) with 6-311++G(d,p) basis set and a SCF convergence criterion of 106 . As the inhibitor molecules are experimentally tested in solution, the optimizations are performed using a polarizable water model (c-PCM) (Barone and Cossi, 1998; Cossi et al., 2002, 2003; Wang and Li, 2009). Further information on the computational details is given in the **Supplementary Material**.

We quantify the structural and chemical similarity between inhibitor structures using the SOAP-REMatch kernel (Bartók et al., 2013; De et al., 2016) to investigate the relation between their structure and associated properties. The SOAP kernel compares local atomic environments and the REMatch (Regularized Entropy Match) kernel condenses the local similarities between two structures into a global similarity measure. A local environment is defined within a spherical region of radius rc centered on an atom and is built by a superposition of Gaussian functions with width ξ . The larger rc is chosen, the more structural information surrounding the atom is included. The SOAP kernel measures the rotationally and translationally invariant overlap between two such local environments and can be raised to a power ζ to discriminate more between large (∼ 0.9) and medium (< 0.6) similarities. The combination of the local similarities can be tuned by the hyper parameter γ of the REMatch kernel. For large values (γ ∼ 10) more equal weights are assigned to the local similarities while for small values (γ ∼ 0.01) only the best matching pairs of local environments are selected to compute the global similarity (see De et al., 2016 for more details).

To help the visualization of potential structure-property relationship, we consider each structure to lie in a highdimensional space defined by the SOAP-REMatch kernel, which is transformed into a distance (Berg et al., 1984), and we project this information on a two-dimensional map using sketch-map (Ceriotti et al., 2011). This dimensionality reduction technique allows to focus the distortions of the space so that close/distant, i.e., similar/dissimilar, structures in the high-dimensional space keep this relationship in the low dimensional space. This behavior is achieved by a sigmoid function that is applied to the distances and is mainly influenced by the switching distance σ, as well as a and b as tuning parameters (see **Supplementary Material**).

Thus, it is possible to create a two-dimensional similarity landscape that allows assessing the molecular similarity by analyzing relative positions and cluster formations. However, due to the form of the sigmoid function, far apart points can be arbitrarily far apart in the lower dimensional projection–making a physical interpretation of distances between basins in the low dimensional projection impossible.

#### 2.3. The Inhibition Prediction Workflow

When combining the presented methods, it is possible to visualize the relationship between the molecular structures and the corresponding inhibition efficiencies in a propertystructure landscape where all experimentally tested structures act as landmark points. Subsequently, the inhibition efficiencies of uninvestigated compounds can be predicted following the proposed workflow (**Figure 1**).

Experimental inhibition efficiencies obtained from a corrosion inhibitor databank (Lamaka et al., 2017), as well as molecular similarity measures determined with the SOAP-REMatch kernel can be combined to create a two-dimensional property-structure landscape for the tested inhibitor molecules. Here, clusters can indicate correlations between the inhibition efficiency and inhibitor structure, allowing to relate certain molecular structures to potentially promoting or inhibiting corrosion properties. For now, the small sample size favors an unsupervised machine learning technique where the decision boundaries are drawn by human, instead of a supervised learning algorithm (Kotsiantis et al., 2007).

Consequently, to predict the inhibition properties of an untested compound its structural relationship to the landmark points has to be determined. This can be accomplished by out-of-sample embedding, where the new structure is projected into the generated sketch-map by reproducing the distances to the previously defined landmark points (Ceriotti et al., 2011). Once the structure is projected into the property-structure landscape, its relative position to previously identified clusters can help to assess the impact on corrosion events. Concurrently, this approach indicates whether it is reasonable to examine untested additives in further corrosion experiments, hence saving a tremendous amount of time and resources compared to an experimental high-throughput approach.

# 3. RESULTS

To create a sketch-map displaying the relationship between corrosion inhibition efficiency and inhibitor structure, a total of 80 compounds was chosen out of the 151 experimentally tested structures provided in a corrosion inhibitor database (Lamaka et al., 2017). All structures were chosen based on a mutual inhibitor concentration of 0.05 M during the hydrogen evolution experiment to avoid concentration dependencies. Before conducting any analysis, the dataset has been further subdivided into 74 plus 6 randomly selected training and test structures – 74 structures for creating the sketch-map and six structures for validating the inhibition workflow.

After geometry optimization, we measure the structural and chemical similarity between these structures using the SOAP-REMatch kernel (De et al., 2016). In order to improve the understanding of the property-structure relationship, the influence of the respective parameters is examined. Indeed, to achieve a wide range of applicability, the SOAP-REMatch kernel and sketch-map technique rely on a few hyper parameters that need to be tuned accordingly (see De et al., 2017; Musil et al., 2018 for a more comprehensive discussion). Depending on the choice of hyper parameters, structural data points are either divided into clusters based on an observable similarity or appear completely scattered. Hence, the parameters have a strong impact on the identifiability of correlations between structure and investigated property. After thorough investigation of the parameter behavior (see **Supplementary Material**), a set of parameters is chosen that allows the division of structural data points with similar corrosion properties into clusters.

To put a higher focus on the local atomic structure for the similarity determination using SOAP, a cutoff radius rc = 3.0 Å is chosen which includes all moieties of interest but neglects the overall molecular structure for most of the investigated dissolution modulators. For a good balance between strict similarity requirements and a sufficient number of pairs of local environments, the Gaussian width is set to ξ = 0.3. By choosing ζ = 2.0 the discrimination between large and medium similarities is increased, thus amplifying clustering effects. Based on the parameter γ = 2.0, a broad selection of well matching pairs of local environments—as determined by the SOAP kernel—is taken into account to compute the global similarity using the SOAP-REMatch kernel.

As the dimensionality reduction with sketch-map is based on a sigmoid function, the corresponding parameters have to be optimized for the given data. Again, optimizing the parameters with respect to cluster formations of the structural data points, choosing σ = 2 as well as the tuning parameters a = b = 3 results in the sketch-maps shown in **Figures 2**, **3**. The data points originating from the input structures are divided into two elongated "islands," a small island in the lower left and a larger island in the upper right part of the sketch-map. It is noteworthy

that aromatic compounds are solely found in the lower region whereas aliphatic compounds are distributed in the upper island of the sketch-map. This leads us to the conclusion, that the chosen parameters are well suited to generate a sketch-map of the investigated molecule database.

When coloring the structural data points according to their corresponding inhibition efficiency (**Figure 2A**), the upper right island is further divided into two clusters, where the left cluster is populated by corrosion inhibitors (green) and the right cluster mostly by corrosion promoters (purple) or moderately inhibiting (light green) additives. The lower left island is dominantly populated by corrosion inhibitors, except for three structures on its outer edge. Cluster formations clearly indicate a propertystructure relationship, allowing to cautiously correlate inhibition efficiency and molecular structure.

For the inhibition prediction it is desired to project not yet tested compounds into the generated sketch-map and relate their position to the three identified clusters. When a new structure is projected into an area with dominantly corrosion inhibitors or promoters, it is assumed to share similar inhibition properties and can be further investigated experimentally if desired. For purposes of validation, six structures of the experimental database are randomly chosen and projected into the sketch-map by determining their global similarity including all structures. Subsequently, the distance of the new structures is related to the 74 defined landmarks and used to compute the required projections. As a guide for the eye, the three previously identified clusters are outlined with dashed lines and colored according to the median inhibition efficiency in the respective region (**Figure 2B**).

As the new structures differ strongly in topology, it is natural that the computed projections lead to differing positions in the sketch-map. In relation to the landmarks, structures containing unusually coordinated atoms, additional atom species or an unusual number of functional groups are projected in regions far away from the observed islands, indicating discrepancies in similarity. However, except for one structure at the top of the sketch-map, structural similarities to the defined landmarks result in projections within or close to the generated sketch-map, where the corresponding inhibition efficiencies match the relative positioning to the inhibitor and promoter clusters fairly well.

The generated sketch-map can also be used to correlate the dissolution modulator structure to other properties, as for instance the HOMO-LUMO gap (HL gap). The energetic difference between the highest occupied and lowest unoccupied molecular orbital (HOMO and LUMO) is indicative for the affinity of the investigated corrosion modulators to transition metals (Griffith and Orgel, 1957), where formation of these complexes is more likely with lower HL gaps as this allows for a energetically more favorable overlap of the involved orbitals. Moreover, the HL gap is a sound indicator for chemical reactivity as the stability of a molecule increases with larger HL gaps. Concomitantly, the reactivity of the dissolution modulator decreases (Aihara, 2000). Consequently, aromatic ligands (e.g., pyridine derivatives) are more likely to form complexes with transition metals (e.g., Fe, Ni) than aliphatic ligands that, in general, exhibit larger HL gaps. Hence, the HL gap might be an important parameter that has to be taken into account in future studies to adequately predict the capability of an untested compound to prevent the re-deposition of noble impurities like

iron (Höche et al., 2016; Lamaka et al., 2016). The HL gaps were calculated on the TPSSh/def2SVP level of density functional theory using Turbomole 7.2 (TURBOMOLE, 2017) for each of the 80 compounds (**Figure 3**). As computing the HL gaps using the B3LYP/6-311++G\*\* level of theory that was employed for the STaGE calculations is computationally rather demanding, TPSSh/def2SVP is chosen here as a fast and accurate alternative. Comparing the optimized geometries for each functional, no structural discrepancies could be observed.

Coloring the sketch-map according to the calculated HL gaps, puts further emphasis on the expected separation between aromatic and aliphatic compounds in the investigated dataset. Aromatic structures in the lower left island are assigned with rather low values of 3.2–5.3 eV whereas aliphatic compounds in the top right island correspond to rather high energy gaps of 5.5–7.4 eV. Albeit this outcome corroborates our current working hypothesis, further work is required to quantitatively correlate the HL gap to the inhibition efficiency of potential inhibitor molecules based on the employed sketch-map approach.

# 4. DISCUSSION

The acquired property-structure landscape in **Figure 2A** uncovers a clear relationship between inhibitor structure and inhibition efficiency, whereas only a few outliers in the defined corrosion inhibitor and promoter clusters are observed. Furthermore, almost all new molecules that are projected into the sketch-map, matching the landmarks in similarity, are correctly positioned within or close to the defined clusters according to their corresponding inhibition efficiency. Hence we are confident, that the presented concept is suitable to predict the potential of uninvestigated corrosion inhibitors or promoters based on their resemblance to a defined landmark structure.