The structure of pomelo peel arouses research interest in recent years because of the outstanding damping and energy dissipating performance of the pomelo peel. Researchers found that pomelo peel has varying pore size through the peel thickness; the pore size gradient is one of the key reasons leading to superior energy dissipation performance of pomelo peel. In this paper, we introduce a method to model pomelo peel bioinspired foams with nonuniform pore distribution. We generate the skeletal open cell structure of the bioinspired foams using Voronoi tessellation. The skeleton of the bioinspired foams is built as three-dimension (3D) beam elements in a full-scale finite element model. The quasi-static and dynamic mechanical behaviors of the pomelo peel bioinspired foams could be derived through a finite element analysis (FEA). We illustrate our method using a case study of pomelo peel bioinspired aluminum foams under quasi-static compression and free fall impact circumstances. The case study results validate our method and demonstrate the superior impact resistance and damping behavior of bioinspired foam with gradient porosity for designers.

## Introduction

The structure of pomelo peel begins to receive research attention in recent years because of the outstanding impact damping and energy dissipating performance of the pomelo peels [1–3]. The pomelo (Citrus maxima), also spelled pummelo or called shaddock, is a fruit native to South and Southeast Asia. A pomelo is able to fall from the height of more than 10 m in its natural environment without cracking or inner damage [1]. Experimental results show that up to 90% impact energy is dissipated during the free fall of a pomelo [4].

The open cell foam structure in a pomelo peel is shown in Fig. 1 [5]. There is a gradient of pore size and shape over the thickness of the peel. The finest pores exist near the outer surface of the peel. The pores become larger and more open along the distance to the outer surface of the peel. In the transition region to inner pulp of the pomelo, the pores change back to fine size but with flat and extended shape. Moreover, the pomelo peel also has many fibrous vascular bundles [2] that act as the struts of the peel. Each bundle includes several interconnected biological cells. The cells are filled with liquids and have multilayer walls.

Designers are encouraged to use a similar varying pore size structure designing open-cell foams [1,6]. However, a mechanical model of bioinspired foam with nonuniform pore distribution that could be used for design has not been published yet. The existing models focused on the open-cell foams with uniformly, periodically, or randomly distributed pores [7–10]. Designers cannot tune the pore distribution to develop a foam material with desired mechanical properties (e.g., impact resistance) through computational simulations.

In this paper, we develop a method to model pomelo peel bioinspired foams with nonuniform pore distribution. We generate the skeletal open cell structure of the bioinspired foams using Voronoi tessellation. The skeleton of the bioinspired foams is built as three-dimension (3D) beam elements in a full-scale finite element model. The quasi-static and dynamic mechanical behaviors of the pomelo peel bioinspired foams could be derived through a finite element analysis (FEA). We illustrate our method using a case study of pomelo peel bioinspired aluminum foams under quasi-static compression and free fall impact circumstances. The quasi-static compression simulation result validates our method since the results we derive for the bioinspired foams have the same general characteristics as those seen in published experimental and theoretical modeling results. The free fall dynamic analysis demonstrates the superior impact resistance and damping behavior of bioinspired foam with gradient porosity.

The case study results show great potential to apply the bioinspired foam with gradient porosity in the areas ranging from simple packaging to design for high speed vehicular impacts. Our method offers designers a useful tool to simulate the mechanical behavior of bioinspired foams under quasi-static or dynamic loading conditions of interest. Design parameters of the bioinspired foam could be optimized in order to obtain desired material properties (e.g., energy dissipation).

## A Pomelo Peel Bioinspired Foam Modeling Method

We introduce a method to model pomelo peel bioinspired foams with nonuniform pore distribution in this section. We generate the skeletal open cell structure of the bioinspired foams using Voronoi tessellation. We program a Python code to read the coordinates of cell vertices (edge points) and remove the duplicate wire segments and the wire segments that are below a tolerance value. The pretreated skeleton of the bioinspired foams is built as 3D beam elements in a full-scale finite element model. The quasi-static and dynamic mechanical behavior of the pomelo peel bioinspired foams could be derived through FEA. To improve the computational efficiency of FEA, we apply a beam to beam contact algorithm and a mass scaling technique.

### Development of Open Cell Structure.

We mimic the pore size gradient and the branching bundles existed in pomelo peels using Voronoi tessellation [11]. An open source C++ library known as Voro++ is used [12]. The Voronoi tessellation includes two steps that are seeding and pore (also called cell) formation as follows.

The algorithm utilizes a Voronoi 3D tessellation of space to create cells of various sizes within a given volume. The program starts by defining the volume that is split into different cells. A set number of points are then randomly placed into the volume, known as seeds. The seeds form the cells of the foam material. Designers can control the seed distribution to construct desired pore size gradient. As stated in Sec. 2, there is more cellular material at the edges of the pomelo peel and less material in the middle of the peel. To replicate the gradient porosity in the pomelo peel, a high number of random seeds can be placed near the top and bottom ends of the volume space followed by a smaller amount of seeds dispersed within the center portion. A similar fashion can be used to other bioinspired application scenarios.

As these seeds are distributed within the foam material, the perpendicular bisector is calculated between one seed and all the other seeds to ultimately form a cell around each seed. The decision-making process is to form a cube around each seed while cutting away portions of it corresponding to the perpendicular bisecting planes calculated with nearby seeds. Eventually, the bisectors lie outside of the cell containing a seed, and the process moves on to another seed.

### Full-Scale Finite Element Modeling.

A full-scale finite element model of the bioinspired foam is established using the coordinates of cell vertices (edge points). We program a Python code to read the coordinates of each point from the Voro++ gnuplot and connect them in sequence if there is no empty line separating the consecutive points to be connected. When we encounter an empty line, this is an indicator that the chain of points is to be broken, and the next point given is not meant to be connected to any previous points, rather only to the points that lie after it (until another empty line is encountered), and the process is repeated until it reaches the last coordinates information.

The duplicate wire segment in the model is then removed through the Python code before we render each line into a 3D beam element. The wire segments that are below a tolerance value also are removed. As the number of wire segments in the model is increased, there are multiple wire segments that have either a length of zero, or some value very close to zero. From the perspective of FEA, if the beams of such short length are permitted in the model, it would cause convergence issues during analysis. To remedy this problem, we remove the segment and merge the segment endpoints if the length of a wire segment is below a given tolerance value. At the same time, we make the appropriate adjustments to the edges connected to those endpoints.

The pretreated wire segment model is then imported into commercial FEA software to construct a full-scale finite element model. Each wire segment is rendered as a 3D beam element in the model. However, importing the wire segments into commercial FEA software as geometric features has shown to take a prohibitive amount of time, because each geometric entity must be rendered and added to the model. To improve efficiency, we write the input file directly from the gnuplot file contents, and the input file is imported into FEA software that produces an orphan mesh. An orphan mesh is not associated with any geometry but instead only has information on the nodal coordinates, element types, and element connectivity. This importation practice has the benefit of greatly reducing the time to import a wire segment model into the FEA software. Once the input file is written, the file is submitted and run without accessing the graphical user interface of the FEA software.

As a last step in importing and conditioning the geometry for FEA, the mesh is queried and any beam elements that are below a given tolerance value are extended until they are at, or above, the tolerance. This treatment would then conclude the conditioning of the geometry and the mesh to ensure reasonable stable time increments and an overall tractable computational problem.

### Quasi-Static and Dynamic Analysis.

The quasi-static or dynamic response of the bioinspired foam under specified boundary condition could be derived using the full-scale finite element model established as above. To improve the computational efficiency of the bioinspired foam FEA, we apply a beam to beam contact algorithm. The contact patch between each two contacting beams is detected and the amount of sliding is computed. The consistent linearization of the contact contribution is computed and added into the governing equations [13]. When using a beam to beam contact algorithm, there is a loss in resolution of the stress occurring at the contact patch between beams, but this is trade-off for a reduction in computational time. The beam to beam contact algorithm is suited for simulation of cellular solids since the simulation is not concerned with high resolution of the stresses at the contact patches between each of the beams. Rather the primary concern of our model is to extract the overall behavior of the foam structure that arises from a large number of contact interactions. The contact algorithm becomes especially useful when introducing many beam elements in the model such that there are thousands of beams exhibiting the nonlinearities associated with contact stresses at any given time [10].

To further improve the computational efficiency of the FEA, we utilize mass scaling on some elements to increase the stable time increment for an explicit FEA [14]. In an explicit analysis, the solution is obtained by integrating through time using many small-time increments that are approximately dictated by the ratio of the element length divided by the dilatational wave speed. Thus, for elements that have a short length compared to others, they are going to have the smallest transit time of the dilatational wave, and therefore, constrain the stable time increment for the entire model. To remedy this problem, the dilatational wave speed is inversely proportional to the density of the material, and thus, the density can be artificially increased in order to decrease the dilatational wave speed and increase the stable time increment. Explicit dynamic procedures involving quasi-static simulations tend to have very long simulation times compared to that of shorter dynamic responses such as impact or shock. Additionally, mass scaling is the preferred method of reducing the solution time if rate dependencies are included in the problem, which is applied in the foam modeling. When used correctly, mass scaling can improve the efficiency of a FEA by allowing the model to run faster while still maintaining a certain degree of accuracy [15,16]. Similar to the import process that remove beams with short length, beams that are above the tolerance but comparatively short tend to constrain the stable time increment of the entire model. Mass scaling can be used uniformly over the entire model or only on elements that are below a target stable time increment. The latter is the method that we suggest for foam modeling. Of note, designers should ensure that the mass scaling does not significantly alter the solution of the FEA, with a primary indicator as a comparison of the internal strain energies to the kinetic energy of the model throughout the analysis.

## Case Study of Pomelo Peel Bioinspired Aluminum Foam

We introduce a method to model bioinspired foam with any given pore distribution in Sec. 2. In this section, we use the quasi-static and dynamic modeling of pomelo peel bioinspired aluminum foam as a case study to illustrate our method. The method could be applied in similar fashion to other application scenarios. The quasi-static compression simulation results validate our method since the results we derive for the bioinspired foams have the same general characteristics as those seen in published experimental and theoretical modeling results. The free fall dynamic analysis demonstrates the superior impact resistance and damping behavior of bioinspired foam with gradient porosity.

### Quasi-Static Compression Analysis.

We use a cylindrical volume to benchmark results with published quasi-static compression experimental data of the pomelo peel. The cylindrical specimen has a height of 3 cm and a diameter of 2.54 cm. We establish two full-scale finite element models with the same seed quantity (8000) following the procedures stated in Sec. 2. The material properties of aluminum 6061-T6 [17] are used since the material data of the pomelo peel bundle are not available. The beams are given a circular cross section with a radius of 0.002 cm according to experimental studies of the pomelo peel [2,3]. The first specimen has uniform porosity as shown in Fig. 2 (left). The second specimen characterizes the gradient porosity of the pomelo peel. We disperse 66% of the pores within 0.6 cm of the top and bottom faces of the specimen as shown in Fig. 2 (right). The average beam length in a specified area is inversely proportional to the seed density, but an individual beam has random length. As stated in Sec. 2.2, any beam that is below 200 *μ*m is extended until the beam is at or above that tolerance.

To replicate an experimental quasi-static compression test, we make the compressive FEA with the two specimens by constraining the bottom face in all directions and rotations and displacing the top face using a rigid body in the shape of a plate. As stated in Sec. 2.3, a beam to beam contact algorithm and mass scaling technique are used to improve the computational efficiency of the FEA.

The quasi-static compressive FEA results of the specimen with uniform porosity are shown in Figs. 3 and 4. The stress–strain curve shows the same general characteristics as those seen in published experimental [18–20] and theoretical modeling [7–10] results of cellular solids with uniform porosity. The three typical regions present in the compaction of cellular materials (elastic, plateau, and densification) can also be seen when studying the stress contours (Fig. 3) during the compaction as well. Initially, there are many beams that are incorporated in resisting the deformation, and then, gradually more beams begin to plastically deform and fail while the beams begin to compact each other.

We also compare the stress–strain curve of the gradient porosity specimen with the curve of the uniform porosity specimen in Fig. 4. We observe that initially both compaction behaviors of the specimens with uniform and changing porosity are very similar, but upon increasing the magnitudes of the reaction forces some differences can be realized. Although the specimen with a changing porosity has the same number of pores and approximately the same number of beam elements, it can produce a smaller reaction force compared to the specimen with uniform porosity.

Figure 4 also shows that our model is able to replicate the stress–strain response of the pomelo peel derived through compression experiments [2,3] although the material properties are different (we use Aluminum 6061-T6). The only location this difference in material properties can be seen is at the initial stages of the loading. There is a notable linear elastic region that coincides with the elastic region of the aluminum material itself. This behavior corresponds to typical compressions of cellular materials in which the initial stages are the beams undergoing simple bending. If the material properties are modified to replicate that of pomelo peel bundles, this elastic regime may disappear. The model would then simulate the behavior of a pomelo peel more accurately.

### Free Fall Impact Analysis.

We also make an impact analysis for the cylindrical specimen with gradient porosity. We simulate real-world conditions subjecting to a free fall from 5 ft (1.52 m) and determine dynamic response of the specimen once the specimen impacts the ground. To simulate the free fall scenario, all the beam elements in the specimen are given a predefined field value of 5.47 m/s. The specimen is then placed 2.1 × 10^{−5} m from a rigid plate that replicate the ground the specimen would impact. The plate is fully fixed in all displacements and rotations. The beam to beam contact algorithm and mass scaling technique are applied in the same way as in Sec. 3.1.

Due to the impacts taking place over very short span of time (5.5 × 10^{−4} s in this case), the analysis for the free fall impact is not as computationally taxing as the prior quasi-static analysis. Figure 5 captures the deformation and stress distributions that take place during the impact time span.

A comparison of the strain energy and the kinetic energy of the bioinspired foam is carried out throughout the duration of the impact simulation to understand how the foam responses to the impact. The comparison is shown in Fig. 6 with an additional energy included which is the plastic dissipation of the foam. It can be seen that the first loading cycle dissipates a high amount of energy, and the second cycles then dissipate much less. This behavior can be attributed to two reasons. The first reason is the plastic deformation of the structure during the loading cycle. Another reason is the frictional dissipation between the beam elements during densification. A notable behavior of the specimen is that the specimen almost immediately enters the densification stage because of the limited free space between the beams.

The impact simulation results above shed light on favorable characteristics of the bioinspired foam with gradient porosity for designers. Figure 6 shows that shock from the impact propagates through the foam but the shock is arrested before reaching the other end. The similar gradient porosity foam design could be used in applications involving high impact, shock, vibration, etc., to protect human or instruments on the other end. Figure 6 also shows that the kinetic energy is turned into strain energy of the structure, while some of the kinetic energy is dissipated through plastic dissipation. After the initial impact stage, the strain energy drops and oscillates around a steady-state value, while part of the strain energy is then converted back to kinetic energy as the specimen at the point it bounces off the ground and is no longer in contact with the ground. The conversion of energies back and forth between strain and kinetic energies ultimately dampen the impact that the foam has with the ground. The energy conversion could serve the purpose of dampening the impact of whatever structure that resides on the other side of the foam specimen.

## Conclusions and Future Work

We introduce a method to model pomelo peel bioinspired foams with nonuniform porosity. We generate the skeletal open cell structure of the bioinspired foam using Voronoi tessellation. The skeleton of the bioinspired foams is built as 3D beam elements in a full-scale finite element model. The quasi-static and dynamic mechanical behaviors of the pomelo peel bioinspired foams can be derived through an FEA. We apply a beam to beam contact algorithm and a mass scaling technique to improve the computational efficiency of FEA. Our method is illustrated and validated through a case study of pomelo peel bioinspired aluminum foams under quasi-static compression and free fall impact circumstances.

The quasi-static compression results show that the bioinspired foam with nonuniform porosity can produce a smaller reaction force compared to the foam with uniform porosity. The free fall dynamic analysis demonstrates the superior impact resistance and damping behavior of bioinspired foam with gradient porosity. These results show that the bioinspired foam could have promising applications ranging from simple packaging to design for high speed vehicular impacts and may also be used for thermal insulation or biomedical purposes. Designers could use our method to tune the design parameters (e.g., pore distribution) in order to derive the bioinspired foam with desired mechanical response under specified boundary conditions through computational simulations.

This paper has several limitations that offer opportunities for future research. We illustrate our method through a case study of pomelo peel bioinspired aluminum foams. More case studies of different foam materials (e.g., nickel and biomaterial) could be carried out by researchers in the future. The foam material with optimized porosity gradient could be fabricated and tested to further validate our method. Moreover, we use Voronoi tessellation to model bioinspired foam with nonuniform porosity. Many other modeling approaches (e.g., splines curves/surfaces, pixel/voxel, periodic surface models, and cellular truss) could be employed and compared with our approach. In addition, the gradient porosity is one of the key features leading to superior energy dissipation performance of pomelo peel. The other features of pomelo peel (e.g., the liquids in the cell) could be explored in future research, and these features may inspire new design concepts.

## Funding Data

National Science Foundation (EFMA-1240483).