Elastic local buckling behaviour of beetle elytron plate

Abstract Beetle elytron plate (BEP) is a biomimetic sandwich structure inspired by the internal architecture of beetle elytra and characterisedby trabeculae in the core. This type of structure has been shown to possess superior mechanical properties to conventional sandwich plates; however, there are no studies that evaluate its structural bending resistance. This paper develops an analytical method to calculate the key component of bending resistance of BEPs: the elastic local buckling load of the compression skin. It assumes that the compression skin of BEPs is simply supported by the trabecular core. After eliminating local buckling in the edges of the compression skin outside the trabeculae, two buckling zones, depending on the ratio ( η ) of trabecular radius to the distance between two adjacent trabeculae, are identified. At low η values ( η ≤ 0.25), elastic buckling occurs in the space of the compression skin surrounded by four adjacent trabeculae. Beyond the critical value of η ( η >0.25), buckling occurs in the compression skin enclosed by individual trabecula. Guided by finite element simulation results, this paper identifies a new suite of deformation shape functions and derives local elastic buckling load for the compression skin according to the principle of minimal total potential energy. Afterwards, a convenient quadratic polynomial regression equation is proposed to modify the elastic buckling coefficient of the compression skin of equivalent conventional grid honeycomb sandwich plates, with the maximum difference between analytical calculation results and finite element simulation results being about 7%.


Introduction
Sandwich plates are a type of structure with high specific strength and stiffness [1][2][3]; therefore, many studies have been conducted to investigate their performance (such as buckling behaviour [4], fire resistance [5], bending and compressive properties [6,7]), energy dissipation and impact resistance [8,9]) and applications in a variety of different industrial sectors (such as aerospace [10], automobiles [11], ships [12] and building structures [13]). As the efficiency of traditional sandwich structures is fully explored, new types of high-performance sandwich structures that feature foam [14,15], truss [16,17] and corrugated cores [18,19] have been developed to satisfy the ever-increasing demand for structural efficiency. Recent studies have also focused on bionic structures, such as nacre [20], biological tissues like fruit walls, muscle, wood and bone [21] and woodpecker beaks [22]. Based upon the latter, Ha et al. [22] proposed a type of sandwich plate with a wavy core with their simulation results suggesting that the new structure has better energy absorbing ability under impact load relative to sandwich structures featuring conventional honeycomb cores.
This research is concerned with a new type of bio-inspired sandwich structure, referred to as the beetle elytron plate (BEP). Over millions of years, beetles have evolved to fly and their body structure is a result  [28]; (b) microstructure of elytron upper skin [28,29]; (c) the trabecular structure [28,29]; (d) biomimetic model, BEP [30]. (Note: broad arrows refer to the trabecular structure; stars refer to the honeycomb wall).
results of these studies confirm that BEP has a higher load-carrying capacity than comparable conventional sandwich plates in bending. However, these studies have so far been limited to phenomenological observations based on experiments. If this innovative structure is to be used in practice, it is necessary to be able to predict its behaviour analytically so that a suitable and rational design methodology can be developed.
The bending resistance of both BEPs and conventional sandwich plates is predominately controlled by the compression resistance of the skin. Since the compression skin is slender, the local buckling load of the skin is expected to be the limiting factor. The first step, which is the focus of this study, is to quantify the elastic local buckling load of the compression skin. Therefore, the objectives of this study are: (1) to understand the different modes of buckling behaviour of the compression skin of BEP, (2) to identify the preferred local buckling mode with the highest buckling load based on simulation work, and (3) to derive an analytical method to calculate the elastic local buckling load. Validated numerical simulations will be used to achieve these objectives.
It should be appreciated that there may be other structural shapes that are more efficient, e.g. using bracing in the core. However, at present, no such solution exists. Therefore, it is not possible to compare the BEP structure of this research with potentially more efficient structures. Nevertheless, more efficient structures will be considered in future research.

BEP structure and research methodology
The basic model of a BEP and its structural parameters are shown in Fig. 2 [30,31]. The BEP consists of two layers of skin with length L, width W and thickness , which are separated by a core structure with a height h. The basic geometrical shape of the core structure can be adjusted to suit functional and manufacturing requirements. The maximum local buckling stress is achieved for a plate aspect ratio of 1 according to [38,39]. The intended construction of BEP is to maximise its resistance. In addition to convenience and without loss of general applicability, a square grid is selected in the present study.
The key characteristic of this biomimetic sandwich plate is the trabecular structure, with radius R, wall thickness and centre-tocentre distance a in both directions, located at the intersections of grid walls. A dimensionless parameter ( = R/a) is used, where = 0 represents a conventional grid plate without trabeculae. The maximum value of is 0.5.
In this paper, the theoretical derivations to calculate the elastic local buckling behaviour of the compressive skin of a BEP are aided by numerical simulations using the general purpose finite element package ABAQUS [40]. Therefore, the research methodology of this paper involves the following main steps: (1) Validation of the general finite element model. There is no direct numerical model for the buckling behaviour of the compressive skin of BEP; the only published experimental studies of the bending behaviour of BEPs [36,37] cannot be used for this validation exercise, because the failure modes were debonding between the skin and the core rather than local buckling of the skin. Therefore, local buckling of the compression skin of the conventional honeycomb sandwich structure is used. Furthermore, the main effect of trabeculae is imposing a circular boundary condition to the rectilinear compression skin. Therefore, elastic buckling of circular plates under compression is also simulated. For validation of ABAQUS simulations, the simulation results are compared against available analytical solutions [38,39].
(2) The validated ABAQUS model is used to perform a parametric study in terms of and b/a ratio to understand the effects of trabecular structure on local buckling behaviour of the compression skin of BEP in Section 4, including local buckling locations and modes, which is then utilised to identify the most effective trabecular structure for maximising local buckling load of the compression skin.
(3) In Section 5, based on the preferred buckling mode, suitable shape functions of the compression skin are selected and used to derive an analytical method based on minimisation of the total potential energy. The results of this theoretical derivation are used to identify the most important parameters which are then used in Section 6 to propose a simple regression equation for application in routine engineering design.

Validation of ABAQUS modelling -comparison against analytical results for conventional sandwich plates
The conventional sandwich plate structure used for validation has the same arrangement as shown in Fig. 2 except there are no trabeculae. As the core is very slender, it is assumed that it can only provide lateral restraint and not rotational restraint to the compression skin. Therefore, the compression skin is modelled as being simply supported (Fig. 3) in this first study of the subject by the authors. Future research will address the effects of rotational stiffness and axial flexibility of the core on the compression skin.
With the above boundary condition to the compression skin, the elastic buckling loads of the exterior region (the green region) and the interior region (the blue region) in Fig. 3 can be calculated using the following equations [38,39]; and the buckling stress of the whole compressive skin should be the lower one from Eqs.
where, refers to buckling stress of the exterior region, is buckling stress of the interior region; m is the number of half waves along the x direction that gives the minimum buckling stress (m = 1 for square interior zones); flexural rigidity per unit width of the plate: The following material properties and dimensions are used: The input material properties (Young's modulus (E) = 5 GPa and Poisson's ratio ((v) = 0.38) are for the purpose of obtaining numerical  results only. These values have been used as they are the values of a type of 3D printing resin [30,32] that may be used in the authors' further research. However, since this study is about deriving a general method to calculate the elastic buckling load, the exact values of E and v do not have any influence on the accuracy of the analytical method.
The thickness of the skin is 1 mm. The core structure is assumed to be square on plan with width and length W = 2a. The edge distance is b. The length of the ABAQUS simulation model is 12a. A wide range of values (10, 20 and 30) for the slenderness parameter = a/t is considered in this study.
According to Eqs. (2) and (3), if b/a<0. 35, local buckling occurs in the interior region of the structures and if b/a>0.35, local buckling occurs in the exterior region. Since these two buckling modes will also occur in BEP structure, therefore, two geometric conditions with b/a = 0.3 and 0.4 are adopted in the numerical simulations.
The structure is modelled using 3-node triangular general-purpose shell (S3) elements with five integration points in the thickness direction. There is no result available for direct validation of numerical simulation of local buckling involving trabeculae. However, the main effect of trabeculae is imposing a circular boundary condition to the compression skin of the BEP. Therefore, additional numerical simulations are carried out for such a compression plate with different slenderness values and the simulation results are compared with the analytical solution in [39], as shown in Fig. 5. The agreement between numerical simulation and analytical solutions can be considered good, with a maximum difference of 4%.

Local buckling behaviour of BEP
Due to the similarity with the conventional grid sandwich structure, only two units of the compressive skin of the BEP structure are presented in Fig. 6. The following dimensions are used: = 1 mm; the ratio of trabecular diameter to trabecular spacing ( = R/a) ranges from 0 to 0.5; the mesh size is 0.5 × 0.5 mm based on the mesh sensitivity study result presented in the previous section.

Buckling modes
There are three potential locations for the occurrence of local buckling in the BEP structure. In addition to the two locations experienced by conventional grid sandwich structures, specifically the space enclosed by the four trabeculae (referred to as interior buckling), and on the edge (termed exterior buckling), it is also possible for local buckling to occur within the circle of the trabecular (named as circular buckling in this paper), as shown in Fig. 6.
In practical design, it is preferable to avoid exterior buckling, because this is an inefficient mode of structural behaviour, which happens when using more materials on the edge would lead to reduced buckling stress. Furthermore, this buckling mode can be easily engineered out by limiting the b/a ratio. For the conventional grid sandwich structure, this value was found to be 0.3, as discussed in the previous section. The existence of trabeculae may change this value. Therefore, an extra parametric study was conducted to investigate the demarking value for BEP structure. Since the thickness of the skin affects the buckling load under these three different buckling modes in the same way, it is expected that changing the skin thickness will not affect the location  of local buckling. Therefore, the parametric study was undertaken for combinations of b/a and R/a ratios. The a/t ratio is taken as 20. Fig. 7 depicts how the location of buckling changes with increasing b/a and R/a ( ) ratios. The buckling locations is determined according to the position of highest amplitude of half wave, such as three examples shown in Fig. 7(a). For the interior and circular buckling, the highest amplitude occurs in the central point of the corresponding region; while it occurs in the middle of the free edge for the exterior buckling. According to Fig. 7(b), it is found that the maximum b/a value that avoids exterior buckling varies slightly with the R/a ( ) ratio.
If the b/a ratio does not exceed 0.3, exterior buckling will not occur. Therefore, in later sections of this paper, the value b/a is set at 0.3.
Similar to eliminating the exterior buckling mode, the end buckling mode [41] can be engineered out, so that there is no compression skin beyond any trabeculae. Therefore, the present study does not consider the end effect.

Effects of trabecular spacing on buckling position
To examine where local buckling occurs (interior or circular), this section will compare ABAQAUS simulation results of buckling stress as a function of R/a ratio for different skin slenderness ( ) values. Due to existence of trabecula in the BEP structure, the following two values may be used to calculate slenderness of the interior buckling area: or: Since Eq. (5) takes into consideration the dimension of trabecular, it is more appropriate and is therefore used in this paper. Fig. 7 presents simulation results of buckling stress as a function of parameter for three slenderness values of 8, 16 and 24. In all cases, the results peak at = 0.3. When < 0.3, elastic local buckling occurs in the area between adjacent trabeculae (Fig. 6, blue area), with increasing buckling load at increasing (see Fig. 8). The maximum local buckling load is about twice of that without trabeculae. However, in the region of = 0.25-0.30, there is little improvement in buckling stress.
When ≥ 0.3, elastic buckling occurs in the compression skin enclosed by individual trabecula (Fig. 6, yellow area), with decreasing buckling load at increasing trabecular radius. Since increasing the   diameter of trabeculae requires more material for the structure, local buckling within the trabecular circle should also be avoided. Therefore, in the remainder of this paper, only buckling in the area bound by four trabeculae will be considered. The conditions are: b/a ≤ 0.3; R/a ≤ 0.25.

Development of the analytical model
The derivation for elastic local buckling load is based on minimisation of the total potential energy.
The strain energy (U ) and the work (V ) done by an axial load ( ) can be calculated by [38,39]: For a given point ( , ) with an infinitesimally small range ± , the equation for calculating the strain energy over this range is: Similarly, the work done by the axial load ( ) becomes: In these equations, where are the actual out-of-plane (z direction) deflections; and is the shape function.
In general, if the structural zone has the same dimensions and same half-wave lengths, unit shape function is used, so the constant ''C'' in Eq. (10) is the maximum deflection of the zone. However, due to existence of trabeculae, the half-wave lengths may be different in different directions and vary at different locations. Therefore, instead of using unit shape function, the half-wave lengths in both directions of any zone are included in the corresponding shape function as detailed below. The value of ''C'' for each zone is then adjusted to ensure that deformation compatibility is maintained in the entire compression skin.
The strain energy and the work done by the axial load of the whole plate are: The total potential energy of a unit compression skin is: For equilibrium, the total potential energy ( ) should be stationary with respect to the degree of freedom. Therefore, the buckling load is

Deformation shapes
In the analytical derivation, only one part bound by four adjacent trabeculae is considered. This is because the same deformation pattern is repeated over all parts, as demonstrated by the numerical simulation result shown in Fig. 7(b).
The results in Fig. 8 indicate that in region of interest for the interior buckling mode ( ≤ 0.25) of compressive skin of BEP, there are three regions. When (=R/a) is low (between 0 and 0.05), the trabeculae have limited effect as the structure changes from the conventional grid structure to BEP structure. This is shown by a very small rate of increase in buckling load as increases. When the value changes in the region between 0.05 and 0.2, rates of increase in local buckling load are larger. The remaining region (between 0.2 and 0.25) has a maximum rate of increase. It means that these three regions of have different deformation behaviours and should be considered separately.
According to [38,39], the local load is determined by the deformation shape and half-wave length. For the compression skin of BEP structure, these two factors are affected by the trabecular structure, as illustrated in Figs. 9-11. For ≤ 0.05, the shape of the whole buckling region can be described by a sine function, which is the same as for the conventional plate, as shown in Fig. 9(a) and Fig. 10(a, b, and c1). For in the range of 0.05< <0.2, near the edges of x (or y) =0 and a/2 (Fig. 10(d)), the buckling deformation can still be described by a sine function, as shown in Fig. 9(b) and Fig. 10(a, b). However, when x or y is equal to R, the deformed shape lies between sine and cosine functions (Fig. 10(c2)). This indicates that at these two positions (Fig. 10(d), dashed boxes) the out-of-plane rotational restraint of the edges has been improved. When is further increased to the range of 0.2 ≤ ≤ 0.25, the corresponding deformation shape can be directly described by a cosine function Fig. 9(c) and (Fig. 10 (c3)).
Such a phenomenon occurs because the arc-shaped edge partially constrains the rotations of the skin in the x and y directions, as shown in Fig. 10(d, dashed boxes). In other words, the rotational restraint on the local edges of the skin is enhanced by the circular geometry of the trabecula. It is also observed that increasing the radius of the arcshaped edge (namely a larger value) enhances the rotational restraint of the skin until the deformed shape resembles a full cosine function.
Therefore, two beneficial effects caused by the trabecular structure can be inferred from the FEM results: (i) a reduced buckling area, and (ii) enhanced rotational restraint. On this basis, the buckling deformation can be categorised into three modes: (i) buckling with a sine function (referred to as the sine buckling mode, Fig. 10(c1)) when ≤ 0.05; (ii) buckling with a cosine function (referred to as the cosine buckling mode, Fig. 10(c3)) when 0.2 ≤ ≤ 0.25; and (iii) buckling mode with a mixture of sine and cosine functions when 0.5< <0.2, Fig. 10(c2). For mode (iii), the exact deformation shape is difficult to be determined accurately. However, the buckling load for this region are expected to be accurately evaluated by assuming the simpler sine deformation shape. Therefore, the exact, but much more complex, deformation shape function was not considered in the present study.
Due to irregularity of the buckling unit in the zone between four trabeculae, it is divided into five zones (I∼V), as shown in Fig. 11. Because circular buckling within each trabecular is prevented at ≤ 0.25, zone V has low deformations, less than 3% that of the maximum buckling deformation value of the buckling area. Furthermore, the ratio of area enclosed in the circle of a trabecular is small compared to the buckling area (less than /16 for ≤ 0.25). Therefore, the energy stored in zone V is ignored for simplicity in the present study.
These observations have informed the detailed proposals described below in Sections 5.2.1 and 5.2.2 for deformation shape functions and half-wave lengths of different regions of the compression zone of BEP structure. The simple lateral supports along the trabecular do not change the in-plane stress flow, therefore, everywhere, the internal linear load acts in the x-direction and is equal to .

≤ 0.05
The deformation shapes of these four zones can be described by sine function. Fig. 12 shows the detailed shape functions of zones I and II. Based on this, the detailed shape functions of the four zones can be determined as follows. Zone I ( ≤ ≤ ∕2 and ≤ ≤ ∕2): where: Zone II ( ≤ ≤ ∕2 and 0 ≤ ≤ ): The deformation shape function is the same as that of zone I in the y-direction. In the x direction, the sine function can still be used, however, for any given point ( , ), due to the existence of the trabecular structure (Fig. 12b), the half-wavelength, , is shortened from a in Eq. (15) to Eq. (18).
Therefore, the shape function in the x direction is: Therefore, for any given point P ( , ), the complete shape function of zone II is: For a given point P ( i , i ), when is equal to R, Eq. (20) for zone II is reduced to Eq. (15) for zone I. Therefore, at the junction between zones I and II, the shape functions (15 and 20) have the same displacement and rotation, thereby satisfying the compatibility condition. Similarly, compatibility is satisfied at other junctions (Eqs. (15) and (21) for zones I and III; Eqs. (20) and (22) for zones II and IV; Eqs. (21) and (22) for zones III and IV). Zone III (0 ≤ ≤ and ≤ ≤ ∕2): Due to the symmetry between Zone II and III, the coordinates of zone II are changed to give the shape function of Zone III: Zone IV (0 ≤ ≤ and 0 ≤ ≤ ): Based on the deformation shapes of zones II and III, the deformation shape function for zone IV is: Fig. 13 shows detailed deformed shapes of zones I and II. Zone I ( ≤ ≤ ∕2 and ≤ ≤ ∕2): In this zone, the half-wave length is a; while the deformed shape changes from sine function on one side to cosine function on the other side. The shape function in the x direction is defined as: in which the coefficients 1 and 2 are defined as: These functions have been selected to satisfy the boundary conditions of: (i) at = , 1 = 1 and 2 = 0; and (ii) at = , 1 = 0 and 2 = 1.
Similarly, in the y-direction: The complete shape function of this zone is: Zone II ( ≤ ≤ ∕2 and 0 ≤ ≤ ): Similar to the transition of deformation shape from one edge to another edge, the shape function along the x direction, 2 (x), is defined  as follows with the consideration in the shortening of the half-wave length.
Because this zone is bounded by a 1/4 circle, the coefficients 1 and 2 are defined as: They are selected to satisfy the boundary conditions at the edges: (i) at = 0, 1 = 1 and 2 = 0; and (ii) at = , 1 = 0 and 2 = 1.
The complete shape function of this zone is: Zone III (0 ≤ ≤ and ≤ ≤ ∕2): For this zone, the coordinates of zone II are swapped to give: in which: Zone IV (0 ≤ ≤ and 0 ≤ ≤ ): For this zone, the shape functions of zone II in the x direction and that of zone III in the y direction are used to give:

Comparison of simplified analytical solution and numerical results
Substituting the general quantities of internal force and deflections in Eqs. (6)- (14) by the specific equations in Eqs. (15) to (38) for deflections will result in a solution of the critical local buckling load ( ) for the compression skin of BEPs. However, the analytical solutions are too complex for application in routine engineering design so a simplified method is sought.
Results from the analytical method reveal that the difference in elastic local buckling loads between the compressive skins of BEP and conventional honeycomb sandwich structure is solely determined by the parameter . Thus, if everything else is the same, it is possible to modify the solution for a traditional grid honeycomb sandwich structure by replacing the constant value of ''4'' in Eq. (39a) by a coefficient k in Eq. (39b) which is a function of . ( By substituting the buckling load of the BEP into Eq. (39), a series of k values is obtained and further fitted as a regress function in terms of .
Fig. 14 compares the modification factor (k) and relationships between the regression equation results and the ABAQUS simulation results for three different compression skin slenderness values, which indicates a very good agreement with the maximum difference being less than 7%.

Conclusions
This paper has presented the results of analytical and numerical studies on the elastic local buckling behaviour of the compressive skin of BEPs. The elastic local buckling load is an important intermediate value for calculating the bending resistance of BEP structures.
The following conclusions can be drawn: (1) There are three possible buckling locations for the compressive skin of BEP structures: interior bucking in the space between four adjacent trabeculae, exterior bucking on the edge between two adjacent trabeculae, and circular buckling within each individual trabecular. When the ratio of the edge distance to trabecular space, b/a, does not exceed 0.3, exterior buckling can be avoided. This is recommended.
(2) The preferred buckling mode is interior buckling, and this happens when the ratio of trabecular radius (R) to trabecular spacing (a) does not exceed 0.25 ( ≤ 0.25). Under this buckling mode, the elastic buckling load of the compression skin of BEPs increases with increasing , reaching a value of about 2 compared to the conventional grid honeycomb sandwich structure, which corresponds to = 0 for the BEP structure.
(3) The deformation shape of the compression skin of BEP structures varies with . This paper has proposed a suite of deformation shape functions for derivation of an analytical solution to calculate the elastic buckling load of the compression skin of BEP structures.
(4) The elastic buckling load of the compression skin of a BEP structure can be obtained by modifying the solution for the compression skin of the equivalent grid honeycomb sandwich structure by replacing the constant value of ''4'' for the latter (Eq. (39)) by a quadratic equation as a function of the parameter (Eq. (40)).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.