The adoption of Nature-inspired strategies to improve materials has fostered the introduction of cavities. But how to mass-produce structures in which a com- plex architecture of cavities is point-to-point fine-tuned to the local and global application requirements? To this aim, we herein report the use of a procedure based on topology optimization and gas foaming. As an example case, a poly- meric foamed beam whose density map is optimized in 3D for three-point bending was designed and produced by gas foaming a purpose-designed preform. The pre- form was produced with polypropylene and by a high pressure autoclave with CO 2as blowing agent. Optical and scanning electron microscopy as well as X-ray microscopy were used to analyze the 3D optimized foamed structures and showed the effectiveness of the foaming design protocol in producing the FEM-optimized structures. A remarkable two-fold increase in the stiffness of the optimized struc- tures was measured with respect to that of the uniformly foamed counterpart 1 with equal overall mass. With the use of a single recyclable material in a sin- gle processing step, this method allows one to conceive the mass production of optimized, therefore more sustainable, plastic parts. Keywords: topology optimization, polymeric foams, gas foaming, structural materials Graphical abstract 2 1 Introduction Bones, tree trunks, and exoskeletons, the par excellence Nature’s structures, are all made by porous materials [1–3] and remodeling is the natural mechanism that lifelong drives the material adaptation to the structural mechanical (and functional) stimuli. As an example, the internal porous architecture of an humerus bone is the optimal and case-specific solution to the complex system of equations describing our walking, running, jumping, etc. [4, 5]. Mimicking the 3D multigraded porous structure (in terms of pore dimension, size, and orientation) of a humerus does not seem difficult today, especially due to modern manufacturing technologies. However the fundamental question is how to (i) design ab initio and then (ii) produce an optimized artificial structure for daily life, such as a car bumper or a serving tray? Are structures in which a complex architecture of cavities is point-to- point fine-tuned to the local and global application requirements conceivable? Topology optimization can be the answer to (i). Formulated by Prager and Rozvany [6] as an alternative to the well-known optimization approaches of size [7] and shape [8], topology optimization first found practical application with the introduction of numerical methods [9] and is now reportedly used in various fields [10, 11]. Topological optimization is often aimed at minimizing the material to be used or maximizing performance (such as stiffness) for a given structural and/or functional application, returning a 3D density map in the reference volume [12]. Similar to natural remodeling processes, the optimization algorithm iteratively removes material where needed, locally reducing body density [13]. In principle, no constraints are imposed on the resulting topology, and therefore the obtained optimal 3D geometry can be arbitrarily complex, both in the external shape and in the internal cavity architecture. Then, producing such an optimal geometry may be a challenge from a manufacturing point of view. Filling the gap between (mathematical) design and production has boosted research fields in design optimization, applied material science, and production tech- nologies [14, 15]. Along this path, the advent of additive manufacturing has played an important role, facilitating the production of complex shape geometries. However, when it comes to production on a large scale, the layer-by-layer process, characteristic of additive manufacturing methods, does not always provide the required production rate. Materials structures and properties can be altered to a large extent by conventional processing, while additive manufacturing only shows limited capability in this instance, often providing parts with diminished properties. Hence, additive manufacturing is only a limited answer to (ii). Polymeric foaming is a group of technologies that is used to massively produce parts, either with simple or complex geometries, possibly through extrusion, injection molding, or steam chest molding. These technologies can be considered environ- mentally friendly and capable of inducing the desired complex 3D geometries and properties. Examples of fields in which polymeric foams dominate the scene, performance- wise and cost-wise, are impact absorption (e.g., in helmets and protective gears), cushioning, and thermal insulation [16–18]. 3 However, foams available on the market and in the scientific and technical literature are spatially uniform in terms of density and pore morphology, apparently with a large space for improvement: graded or multigraded foams, Nature tells us, could further improve actual performances toward more efficient, optimized use of the matter. Recently, a technology has been introduced to produce multigrade polymeric foams based on time-varying boundary conditions in gas foaming [19]. It has been exploited to produce samples with simple geometries, such as slabs or beads, and has proven to be efficient in targeting specific compressive properties of sintered foamed beads [20]. Accordingly, in the present work we introduce a novel procedure to design and produce a topologically optimized 3D graded-foam structure through physical foaming. With the knowledge of the load conditions and constraints, as well as the mechanical properties of the expanded polymer, an optimized 3D density map is calculated with a topology optimization module within a FEM software. The density map is then produced by gas foaming from a suitably designed preform. Optical and scanning electron microscopy as well as X-ray microscopy (XRM) are used to analyze newly developed foams to verify the effectiveness of the production method. The 3D optimized foam structure achieved is also mechanically validated and compared with spatially uniform foams with the same overall mass. In a single processing step, this novel methodology allows one to conceive the mass production of optimized plastic parts with a single recyclable material, reducing both costs and material usage for a sustainable use of the matter. 2 Results Following the workflow depicted in Figure 1, we first designed an optimized 3D den- sity map ( Design section), which has been then produced via gas-foaming ( Production section) and mechanically tested in the three-point bending (TPB) configuration (Validation section). 2.1 Design, computational topology optimization The result of the computational topology optimization procedure gives an optimized 3D density map, as shown in Figure 2 (see subsection 4.2.1). The density ranges from 100 to 900 kg/m3, lower bounded by the foaming yield of the specific adopted polymer (polypropylene, PP) and upper bounded by the density of the neat, unfoamed polymer. As expected, based on the load conditions and the known results of structural mechanics, the virtual optimized 3D density map reveals a bridge-like shape of the dense portions of the beam ( ∼900 kg/m3, red-colored in the figure). Tinted blue indi- cates regions of lower density that experience minimal stress. The transition from red to blue is gradual, as required in view of the selected physical penalization coefficient defined in the context of the adopted topology optimization approach (see subsection 4.2.1). 4 Fig. 1 The workflow, divided in three phases: design, production and testing. In the design phase, the load configuration, constraint and domain are set and a virtual optimized 3D density map is obtained by the use of a topology optimization computational tool. In the production phase, a preform is conceived and gas foamed in a box-shaped mold. In the testing phase, the obtained foamed beam endowing the physical optimized 3D density map is investigated either in terms of morphological and of structural features. 2.2 Production, optimized foamed beams Gas foaming technology can be used to convert the virtual optimized 3D density map into a physical 3D structure. The preliminary production of a preform is necessary, which is subsequently foamed by the high-pressure quench method. Figure 1 shows the image of the preform, constructed from PP slabs of varying thicknesses, which have been formed into a specific shape (see subsection 4.2.2). By adopting a non-equilibrium gas sorption step of the gas foaming procedure, the gas concentration in the polymer prior to foaming is set pointwise in the preform. At the pressure quench, the polymer foams to the extent determined by the local gas concentration according to the virtual optimized 3D density map. Figure 1 reports the image of a foamed beam. 5 Fig. 2 Virtual optimized 3D density map by computational topology optimization. On the left, the axonometric projection and top, front, and bottom views. On the right, map slices at different distances from the center, where contour plot is used to emphasize isodensity lines. 2.3 Validation The efficacy of the proposed procedure to obtain foamed beams with the optimized 3D density map is assessed by morphological and structural analysis. 2.3.1 Morphology Results of 2D (by scanning electron microscopy, SEM) and 3D (XRM) imaging of the foamed beams are reported in Figure 3 in comparison with the virtual optimized 3D density map. In detail, Figure 3 (a) reports selected sections of the XRM analy- sis corresponding to the images reported in Figure 2. Here, the bright areas represent the polymer, and the black areas correspond to the voids. The microscopic resolu- tion of the XRM gives details of the porous structure. A volume-average procedure schematized by means of Figure 3 (c) (see subsection 4.2.6) returns a density map such as the one reported in Figure 3 (d) in comparison with the corresponding virtual density map of Figure 3 (b). The physical optimized 3D density map resembles the virtual one, in particular in bringing the dense, structural parts in the desired position, proving the effectiveness of the proposed production method. The manual preform production represents an obvious improving direction. Alternative, high-precision, and high-throughput technologies like injection molding could be adopted. The SEM images of Figure 3 (e) reveal details of the density and morphology of the graded foam. 6 Fig. 3 (a) XRM scans at different distances from the center. (b) The 3D virtual optimized density map section at 46.5 mm (from Figure 2). (c) enlarged image from the XRM of the section at 46.5 mm with areas selected for SEM analysis. (d) results of the volume-averaging procedure on the same image. (e) graded foam details by SEM at different magnification. 2.3.2 Three-point bending test Structural properties of foamed beams are measured with TPB tests (Figure 4(a), see subsection 4.2.7 for details): topology optimized and spatially uniform PP foamed beams are compared to prove the value of the proposed approach. Uniform foams were produced from rectangular slabs as detailed in subsection 4.2.4. The load-displacement curves are shown in Figure 4 (b), which shows markedly superior properties of the topologically optimized foams with respect to the uniform ones. The flexural stiffness calculated from the above curves is shown in Figure 4(c) (refer to subsection 4.2.8). The flexural stiffness of spatially uniform foams with a density of 225 kg/m3is 3.1±0.5 N/mm, while for topology optimized foamed beams (with equal mass, i.e., equal apparent density), the flexural stiffness is 6 .5±0.7 N/mm, with a remarkable 210% increase. The optimization procedure is also carried out to produce reduced mass, topologically optimized foams, which gave a flexural stiffness of 4.0±0.5 N/mm in the case of 82% of mass and a flexural stiffness of 1 .9±0.3 N/mm for the 64% of mass. As demonstrated in Figure 4(d), uniform foams can be replaced by topologically optimized ones, resulting in a 26% reduction in mass while maintaining the same performance. These results confirm that the proposed approach is suitable for producing topo- logically optimized foams via physical foaming with enhanced mechanical properties with respect to uniform foams when load conditions are known. 7 Fig. 4 (a) Experimental setup for TPB test. (b) mechanical behavior in TPB configuration for spatially-uniform and topologically optimized foams with equal overall mass. Shaded areas represent the full range of experimental load-displacement curves, while continuous lines represent the average behavior. (c) boxplot and statistical distribution details of stiffness values for topologically optimized foams and spatially-uniform foams. (d) stiffness values as a function of mass percentage with respect to the mass of the uniform foam. Dashed the linear best-fit. 3 Conclusions Industries are now facing the issue of mass-producing structures with a complex archi- tecture to meet specific functional and/or structural application requirements and yet reduce material consumption. Here, we propose a foaming method to efficiently produce optimized polymeric foam structures with any complex architecture of the density and porosity. Taking as a loading condition example the TPB test and using topology opti- mization software, we obtain the optimized 3D density map. From the optimized 3D density map, we produce a preform that, once foamed, gives a foamed sample with the desired 3D density map. Good matching is achieved between the virtual and the physical optimized 3D density map, as confirmed by XRM and SEM observations. Mechanical testing reveals the efficacy of the proposed method, with a two-fold increase in the flexural stiffness of the topology-optimized foamed beams relative to the uniformly foamed beam with equal mass. 8 4 Materials and methods 4.1 Materials Polypropylene (PP) with a density of 900 kg/m3is supplied by the Sinopec Zhenhai Refning & Chemical Company (Zhejiang Province, China) in the form of beads. CO 2 is provided by SOL S.p.A. (Monza, Italy). Compression molded PP slabs were shaped and foamed to produce foamed samples with different geometry for respective testing. 4.2 Methods 4.2.1 Topology optimization Topology optimization is carried out using the Solid Mechanics and Optimization Module based on a Solid Isotropic Material with Penalization (SIMP) approach [9, 21], within COMSOL Multiphysics software. The SIMP approach introduces a continuous relative density control variable, θc(x)∈[θmin,1]⊂R, with x∈Ω⊂R3, where Ω represents the object domain, in which loads and boundary conditions are defined; in particular, boundary tensions, t, act on ∂Ωt(i.e., on a subset of the boundary of Ω), while body forces, b, act on Ω. Furthermore, the SIMP approach introduces a lower limit condition on the variable θc(x), i.e., θc(x)≥θmin>0, and this requirement is introduced to guarantee FEM stability when solving for equilibrium [22]. The relation between the material Young’s modulus, ESIMP (x), and θc(x) is then expressed by ESIMP (x) =Es[θmin+ (1−θmin)θp c(x)], (1) where Esis the Young’s modulus of the dense material and pis a penalization coef- ficient. As classically done in the SIMP approach, we set the Poisson ratio ν0to be independent of θc. In literature, the penalization coefficient pin equation (1) is chosen to properly penalize Young’s moduli, without incurring in numerical difficulties; a classical choice isp= 3. On the basis of experimental evidence, we use p= 2.52, see Appendix A. Topology optimization is driven by an optimization criterion. In the following, we choose the minimization of compliance, while imposing a constraint on the total mass to be distributed within the domain Ω. Accordingly, assuming a linear elastic material 9 behavior, the problem to be solved becomes   minimizeZ Ωb·udΩ +Z ∂Ωt·uds subject toZ ΩC(θc(x)) :ϵ(u) :ϵ(v)dΩ =Z Ωb·vdΩ +Z ∂Ωt·vds subject toZ Ωθc(x)dΩ Z ΩdΩ=Vfrac(2) with the equilibrium written in a weak variational form through equation (2.2) and where C(θc(x)) is the 4th-order stiffness tensor of the elastic material, ϵis the elastic strain tensor, Vfracis the target total volume fraction of the material to be distributed within the domain Ω, uis the displacement and vis a displacement test function. The software COMSOL Multiphysics solves the minimum problem (2) iteratively by adjusting the control variable θc(x) distribution [23]. The solution algorithm measures convergence in terms of the total elastic energy WE, defined as WE=1 2Z ΩC(θc(x)) :ϵ(u) :ϵ(u)dΩ. The spatial domain object of the optimization procedure is selected based on ASTM D790-17 prescription [24–27] on TPB beam shape and on the capacity of adopted foaming equipment. In our case, dimensions are l= 93 mm, b= 19.3 mm and d= 7 mm. The span length Lis 80 mm, as shown in Figure 5 (a). 4.2.2 Preform production To obtain preforms, both for topologically optimized and spatially-uniform foams, PP slabs are used as building blocks. Pieces are cut from PP slabs produced starting from PP beads by using a compression molding machine (model P300P, Collin Gmbh, Ebersberg, Germany). The slabs were produced as follows. The mold is pre-heated to 200 °C, pellets are then placed on the flat press table for 2 min without imposing any pressure. Then, pressure increased up to 5 bar in 1 min, then up to 10 bar in 1 min, and up to 20 bar in 1 min. After the plates cooling (to 25 °C), we extract the PP sheet. The preform for the uniform foam is a slab 2 mm thick with a weight of 2 .8 g. After foaming, the average density is equal to 225 kg/m3. The topologically optimized preform is then designed to obtain, after foaming, an optimized foamed beam with the optimized 3D density map according to the topology optimization result. First, the optimized 3D density map is divided into three regions, 10 each one corresponding to an interval of densities: ρ= 850 −900 kg/m3(first region), ρ= 300 −850 kg/m3(second region) and ρ= 100 −300 kg/m3(third region). For partially/fully foam slabs, the gas penetration depth is determined knowing the PP-CO 2mutual diffusivity, DPP-CO 2(in the hypothesis of Fickian diffusion), that is of the order of 10−9m2/s [28–30]. Fixing the saturation time for the optimized preform equal to ts= 120 s (see subsection 4.2.4), the gas penetration depth is xCO2∼p DPP-CO 2ts∼10−4m. To design the preform, the first region is substituted by the bulk material (i.e., the density of the first region is approximated to be 900 kg/m3), with 2 mm thick slabs. In this case, during the expansion phase, the material remains mostly unfoamed. In the second region, the region that contains the density gradient, PP slabs have a thickness of 1 mm, in order to obtain a graded density distribution after the pressure release. In the third region, slabs with a thickness of 0 .4 mm are used: the density will reaches its lowest value for the adopted foaming conditions at the expansion step due to a complete saturation of slabs and a large volume available in the mold. Three types of optimized preform are produced: with a weight of 2 .8 g (i.e., the same mass of the uniform preform), of 2 .3 g (82% of the mass with respect to the uniform preform) and 1 .7 g (61% of mass with respect to the uniform preform). After the expansion step, the apparent average density of the 2 .8 g optimized foam is equal to 225 kg/m3. Figure 5 (b) reports the image of the box-shape mold, the topologically optimized (2.8 g) and spatially-uniform preforms. 4.2.3 Mold filling validation Mold filling plays a crucial role in determining the final morphology of the foam, and hence the resulting mechanical properties. Foam production prediction of mold filling is needed. The software Moldex3D®is used for such a purpose. The simulation is carried out using the injection molding module. To implement the foam flow front in time, a preliminary foaming experiment is performed starting with a single slab (Figure 5 (c) upper part). A good approximation for volume expansion, V(t), is given by: V(t) =L(t)Al (3) where L(t) is the flow front length in time and Alis the slab surface. L(t) is measured from the experiment via image analysis with the software ImageJ®. The resulting flow rate, Q(t), evaluated with the equation Q(t) =d dtV(t) (4) is implemented in Moldex3D®and controls the optimized preform expansion (Figure 5 (c) lower part) in the mold. 11 4.2.4 Gas foaming The batch foaming equipment is a 0 .3 L custom-made thermo-regulated pressure vessel (model BC-1, High Pressure Equipment Co., Erie, PA, USA), as described in detail elsewhere [31]. The electrical heater is controlled by a PID thermoregulator (mod. 1850, Gefran S.p.A., Provaglio d’Iseo, Brescia, Italy), which takes as input the temperature inside the vessel and the temperature of the heating jacket by using a Pt100 and a J-thermocouple, respectively. A pressure transducer (TK, Gefran S.p.A., Provaglio d’Iseo, Brescia, Italy) is used to measure pressure during the saturation step. The pressure discharge system consists of a discharge ball valve (model 15-71 NFB, High Pressure Equipment Co.,Erie, PA, USA) and an electromechanical actuator (model 15-72 NFB TSR 8, High Pressure Equipment Co., Erie, PA,USA). The pressure drop rate (PDR) is measured through a PLC (Siemens, SIMATIC S7-1200 CPU) with a sampling rate τm= 1 ms. To obtain foamed samples, PP preforms (both for uniform and for topologically optimized) were placed in a rectangular mold with the same dimensions of the reference domain of the topology optimization simulation. The mold with the preform was placed inside the pressure vessel, and the melt temperature, 156 °C, was reached under vacuum conditions. Saturation was performed at 144 .5°C and 120 bar for 300 s in case of preform for uniform foam. In the case of topology optimized preform, saturation was performed at the same pressure and temperature for 120 s. Foaming temperature was 120 °C for both preforms. After the foaming process of both types of preforms, produced foams achieved the same apparent average density equal to 225 kg/m3. 4.2.5 Density measurements Density of TPB samples is measured according to ASTM D792-20 with an analytical balance (MS Semi-Micro Model, Mettler Toledo, Milan, Italy). 4.2.6 Morphology characterization Microstructures of foams are observed either by scanning electron microscopy (SEM) for 2D image analysis and by X-ray microscopy (XRM) for the 3D view. In the case of SEM, foam sections are coated with gold (208HR high-resolution sputter coater, Cressington Scientific Instruments, Watford, England) before imaging with a field-emission-gun-scanning electron microscope (FEG-SEM; Carl Zeiss AG, Oberkochen, Germany). The overall macroscopic 3D-view of the foam microstructure and morphology is obtained via XRM (SKYSCAN 1272, Bruker, Belgium). Experimental density maps are obtained through a volume-averaging procedure starting from XRM scans by using MATLAB®. The volume-averaging procedure consists of number averaging the inten- sity (on a 0-255 scale) of 75 ×75×75 pixels (corresponding to 0 .1 mm3) and normalizing to 900 kg/m3. 12 4.2.7 Three-point bending test The test is performed using an Electromechanical Universal Testing Machine (MTS Insight 10, MTS System Corporation) with a 250 N load cell. The tests are carried out at room temperature (22 ±2°C) in a displacement control setting at a displacement rate of 3 mm/min. The span length is L= 80 mm. The overhang (6 .5 mm per side) is experimentally proved to be sufficient to prevent the specimen from sliding from the supports. 4.2.8 Young’s Modulus and stiffness derivation TPB and compression tests (Appendix A) show that PP foams have a strongly non- linear response, therefore the linear part of both types of curves is determined as follows. The only hypothesis is that there exists a region in which the material behave as linear in the load-displacement curve, and is therefore of the form F=kd, (5) where Fis the force and dis the displacement. If present, the TOE region (i.e., sample adjustment between plates) is eliminated before the data analysis and the load-displacement curve is re-scaled ( F(d) = 0 ↔ d= 0). The load-displacement curve is then fitted in a least-square sense with MATLAB® by using Equation 5 increasing step-by-step the number of displacement points used for the fitting. At the same time, the average error on the measure data ( eF) is computed as follows for each number of acquired points ( N): eF((Fexp1, dexp1)...,(Fexpi, dexpi), ...(FexpN, dexpN), kfit) =... ...=vuut1 N−2NX i=1 Fexpi−kfitdexpi2,(6) where kfitis the fitting parameter (and it depends on N),dexpiis the i-th imposed deformation and Fexpiis the corresponding measured force. A FexpN−eF(FexpN) curve is obtained. Given FexpN, this curve should represent the error associated to the approximation of the material behavior as linear until FexpNis reached. If the behavior is truly linear, it is expected that this error is comparable with the error on the measurement given by the load cell, eLC(FexpN), which is available from the load cell calibration data sheet. In other words, eFis assumed to be a measure of the error of the linear behavior hypothesis. Therefore, we assume that in the region for which eF(FexpN)< eLC(FexpN) the material behavior is linear, while it is non-linear outside. This implies that the non-linearity starts when eF(F∗ expN) =eLC(F∗ expN). In the case of TPB load-displacement curves, the stiffness is obtained by fitting the curve with Equation 5 before F∗ expN. 13 In the case of compression load-displacement curves, the Young’s modulus, Ef, is obtained in the same way but imposing Ef=kh/A , where handAare the height and the area of the compression sample, respectively. Fig. 5 (a) Topology optimization domain: 3D view, top view and front view. (b) Mold for the gas foaming process, along with the PP slab (i.e., the preform to get a uniform foam) and the PP optimized preform (i.e., the preform to get the spatially-optimized graded foam). (c) Mold filling modeling. We used a pilot experiment (upper part) to determine the expansion rate of a slab during foaming to calibrate model parameters. Then, we simulate the preform foaming (lower part). 4.2.9 Data analysis Data analysis is performed by using OriginPro®and MATLAB®. When the same sub- set of measures is considered (for repeatability), data are subjected to the Chauvenet criterion and outliers are rejected. 14 Appendix A Experimental evaluation of the penalization coefficient, p In the standard literature, the penalization coefficient p, introduced in subsection 4.2.1 within the topology optimization procedure to properly penalizing the Young’s modulus when density decreases, is chosen in such a way to avoid numerical difficulties when solving the global minimization problem. A classical choice is p= 3. It is interesting to highlight that Lipka and Ramm [32] noted that the relation gen- erally proposed in the SIMP approach between ESIMP andθcbased only on numerical stability is very similar to the experimental one established between the Young’s mod- ulus of an open-cell foam, Ef, and its density, ρ, as proposed by Gibson and Ashby [33]: Ef=CEsρ ρsn (A1) where ρsis the density of the full-dense (solid) material, Cis a constant (experi- mentally evaluated as C≈1) and nis the so-called Ashby exponent (experimentally evaluated in the range between 1 .5 and 2 .5). Following relation (A1), we measured nfor the material adopted herein, PP, by compression testing several foams with different relative densities. We emphasize that imposing p=n, the relative density map obtained by simula- tionθc(x)|mincoincides with the relative density map ρ(x)/ρs, i.e., the density that should be spatially assigned to a foam to reproduce the optimized ESIMP (x) map. To obtain foamed samples with different densities for the compression test, several solid PP disks with different geometries are used as preforms and placed in cylindrical molds of 9 mm height and 23 mm diameter. Samples are produced with the foaming procedure reported in subsection 4.2.4 and a sorption time of 420 s. Compression tests are performed using an Electromechanical Universal Testing Machine (MTS Insight 10, MTS System Corporation) with a 10 kN load cell. Tests are carried out in displacement control with a displacement rate of 3 mm/min. Ten different groups of densities are produced and tested in three replicates. Den- sity is measured according to ASTM D1622-03 with the same analytical balance of subsection, and the Young’s Moduli are obtained as described in subsection 4.2.8). Assuming C= 1, consistently with equation (1) and experimental observation, the fitting procedure gives n= 2.52±0.04 (see Figure A1). 15 Fig. A1 Derivation of the Ashby exponent n. (i) PP disk-preforms with the same thickness and different (increasing) diameters; (ii) foamed samples with different (increasing) densities produced via physical gas foaming and tested in compression to compute Young’s moduli E; (iii) Normalized Young’s modulus, E/E s, versus density ratio, ρ/ρsin a bi-logarithmic scale, and linear-fitting. Acknowledgments. The X-ray microscopy (XRM) scans and analysis are kindly provided by Dr. Mostafa Borji, Bruker Belgium (Kontich, Belgium). Data Availability The data that support the findings of this study are openly available in figshare, DOI: 10.6084/m9.figshare.24119583.