Multi-objective structural optimization of honeycomb cells

1. Introduction

The metallic and polymeric honeycomb cores widely used in panel structures, energy absorbers and insulation layers are formed by periodic permutation and combination of prismatic cellular microstructures [1-3]. Their tensile/compression and shear properties can be improved effectively by the cell structure optimization [4, 5]. Hence, much effort is usually spent on the optimization of geometric parameters of the cell.

Qin et al. [6] optimized the cell wall thickness and its variation gradient index using Non-dominated Sorting Genetic Algorithm (NSGA) to maximize the energy absorption ratio and minimize the collision force peak. Zhang et al. [7] optimized the cell size of a honeycomb structure by NSGA to maximize the energy absorption efficiency. Qiu et al. [8] optimized the size and key point coordinates of a hexagonal cell using Genetic Algorithm (GA) to minimize the flexibility of a chiral honeycomb panel. Wang et al. [9] optimized the relative density and size of the cell to improve the heat dissipation capacity of a composite sandwich honeycomb panel. Zhang et al. [10] optimized the cell wall thickness of an aluminum honeycomb using software Hyperstudy and LS-DYNA integrated optimization technique. Sorohan et al. [11] adopted the optimization module of ANSYS software to optimize the cell geometric parameters of a hexagonal honeycomb core with out-of-plane isotropic properties. Cinar et al. [12] optimized the thickness, diameter and inner angle of a three-dimensional cell using GA combined with a finite element solver to minimize the natural frequency difference square sum of the honeycomb panel. Gholami et al. [13] and Namvar et al. [14] optimized the cell structure parameters by Particle Swarm Optimization and GA to minimize the deflection of a honeycomb core under uniformly distributed normal load and static load, respectively.

To date, it has not been investigated that on the premise of ensuring the honeycomb core light weighting, the cell structure is optimized to improve the tensile/compression and shear moduli in in-plane. The present study firstly derives the formulas of the equivalent elastic moduli of three types of cells by Energy Method, and analyzes the calculation accuracy by comparison with the results of Refs. [15] and [16]. Then, the multi-objective optimization model is established to maximize the equivalent elastic moduli, and the size and shape of the cell are optimized using GA. Finally, the structural performance of original and optimized honeycomb cores are investigated using statics and modal analysis.

2. Equivalent elastic modulus of the cell

The inclined wall length $l$, vertical wall length $\alpha l$, wall thickness $\beta l$ and inner angle $\theta$ of the basic cell, and the honeycomb core thickness $\eta l$ are shown in Fig. 1. The non-dimensional treatment for the inclined wall length is produced, that is, $\left[l\right]=1$, and the dimensionless quantities of the vertical wall length, wall thickness and honeycomb core thickness are $\alpha$, $\beta$, $\eta$, respectively.

According to $\theta >0°$, $\theta =0°$ and $\theta <0°$, the cell is divided into three types: hexagon cell, quadrilateral cell and concave hexagon cell, and the quarter structure and size parameters are shown in Fig. 2. The tensile and compressive stress in the $X$ or $Y$ direction and shear stress in the $XY$ plane generate when the force acts in the honeycomb core plane. The equivalent elastic and shear moduli are derived, and the processes are the same for hexagon, quadrilateral and concave hexagon cells. The concave hexagonal cell is chosen as a testbed to illustrate the calculation of material constants.

Fig. 1Honeycomb core

a) Honeycomb core

b) Basic cell

Fig. 2Quarter structure and size parameters of hexagonal, quadrilateral and concave hexagonal cells

a)$\theta >0°$

b)$\theta =0°$

c)$\theta <0°$

2.1. Equivalent elastic modulus in the $\mathbf{X}$ direction

The basic cell is subjected to a normal stress along the $X$ direction, ${\sigma }_X$, when the tensile deformation of the honeycomb core occurs, as shown in Fig. 3(a). The displacement of a quarter of the cell is calculated due to the symmetry, and is the same as that of the inclined wall under the combined action of stretching and bending. The vertical wall is regarded as a fixed-end constraint in view of the interaction between cells, and then the inclined wall is simplified to be a cantilever beam model, as shown in Fig. 3(b). The tensile load $F_X$ and bending load $M_X$ are expressed by:

1

$\left\{\begin{array}{l}{F}_X=\frac{1}{2}{\sigma }_X\alpha \eta l^2,\\ M_X=\frac{1}{4}{\sigma }_X\alpha \eta l^3\mathrm{s}\mathrm{i}\mathrm{n}\theta .\end{array}\right.$

Fig. 3Concave hexagonal cell subjected to σX and its equivalent model.

a) Concave hexagonal cell subjected to ${\sigma }_X$

b) Equivalent cantilever beam model

The axial force $F_{XN}\left(x\right)$ and bending moment $M_{XZ}\left(x\right)$ of the inclined wall OA are as follows:

2

$\left\{\begin{array}{l}{F}_{X\mathrm{N}}\left(x\right)=\frac{1}{2}{\sigma }_X\alpha \eta l^2\mathrm{c}\mathrm{o}\mathrm{s}\theta ,\\ M_{XZ}\left(x\right)=\left(\frac{1}{4}l-\frac{x}{2}\right){\sigma }_X\alpha \eta l^2\mathrm{s}\mathrm{i}\mathrm{n}\theta .\end{array}\right.$

The total deformation energy under the coupling action of $F_{X\mathrm{N}}\left(x\right)$ and $M_{XZ}\left(x\right)$, $U_X$, is as follows:

3

$U_X=U_{M_{{\mathrm{}}_{XZ}}}+U_{F_{\mathrm{N}}}={\int }_0^l\frac{{M_{XZ}}^2\left(x\right)}{2E_{\mathrm{s}}{I}_Z}dx+{\int }_0^l\frac{{F_{X\mathrm{N}}}^2\left(x\right)}{2E_{\mathrm{s}}A}dx=\frac{{{\sigma }_X}^2{\eta }^2{l}^7\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \mathrm{c}\mathrm{o}{\mathrm{s}}^4\theta }{32E_{\mathrm{s}}{I}_Z}+\frac{{{\sigma }_X}^2{\eta }^2{l}^5\mathrm{s}\mathrm{i}{\mathrm{n}}^4\theta \mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta }{2E_{\mathrm{s}}A},$

where $U_{M_{XZ}}$and $U_{F_{XN}}$ are the deformation energies produced by $M_{XZ}\left(x\right)$ and $F_{XN}\left(x\right)$, respectively; $E_s$ is the elastic modulus of matrix material; $I_Z$ and $A$ are the moment of inertia and the area of the section.

The displacement of the cell along the $X$ direction, ${\delta }_X$, is obtained by Castigliano’s Theorem:

4

${\delta }_X=\frac{\partial U_X}{\partial F_X}=\frac{{\sigma }_X\eta l^5\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \mathrm{c}\mathrm{o}\mathrm{s}\theta }{12E_s{I}_Z}+\frac{{\sigma }_X\eta l^3\mathrm{c}\mathrm{o}{\mathrm{s}}^3\theta }{E_{s}A}.$

The strain along the $X$ direction, ${\epsilon }_X$, is as follows:

5

${\epsilon }_X=\frac{{\delta }_X}{\alpha l}=\frac{2{\sigma }_X\mathrm{c}\mathrm{o}\mathrm{s}\theta (\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta +{\beta }^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta )}{\alpha {\beta }^3{E}_s}.$

The equivalent elastic modulus along the $X$ direction, $E_X^{*}$, is expressed according to Hooke’s Law:

6

$E_X^{*}=\frac{{\sigma }_X}{{\epsilon }_X}=E_s\frac{\alpha {\beta }^3}{2\mathrm{c}\mathrm{o}\mathrm{s}\theta (\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta +{\beta }^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta )},\left(\theta <0°\right).$

For hexagonal and quadrilateral cells, the equivalent elastic modulus along the $X$ direction is written as:

7

$E_X^{*}=E_s\frac{(\alpha +2\mathrm{c}\mathrm{o}\mathrm{s}\theta ){\beta }^3}{2\mathrm{c}\mathrm{o}\mathrm{s}\theta (\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta +{\beta }^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta )},\left(\theta \ge 0°\right).$

2.2. Equivalent elastic modulus in the Y direction

The equivalent elastic modulus in the $Y$ direction is calculated based on the cantilever beam model when the cell is subjected to a normal stress along the $Y$ direction, ${\sigma }_Y$, as shown in Fig. 4. The tensile load $F_Y$ and bending load $M_Z$ are expressed by:

8

$\left\{\begin{array}{l}{F}_Y={\sigma }_Y\eta l^2\mathrm{c}\mathrm{o}\mathrm{s}\theta ,\\ M_Z=\frac{1}{2}{\sigma }_Y\alpha \eta l^3\mathrm{c}\mathrm{o}\mathrm{s}\theta .\end{array}\right.$

Fig. 4Concave hexagonal cell subjected to σY and its equivalent model

a) Concave hexagonal cell subjected to ${\sigma }_Y$

b) Equivalent cantilever beam model

The axial force $F_{YN}\left(x\right)$ and bending moment $M_{YZ}$ of the inclined wall OA are as follows:

9

$\left\{\begin{array}{l}{F}_{YN}\left(x\right)={\sigma }_Y\eta l^2\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\theta ,\\ M_{YZ}\left(x\right)=\left(\frac{1}{2}l-x\right){\sigma }_Y\eta l^2\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta .\end{array}\right.$

The total deformation energy under the coupling action of $F_{YN}\left(x\right)$ and $M_{YZ}$, $U_Y$, is as follows:

10

$U_Y=U_{M_{YZ}}+U_{F_N}={\int }_0^l\frac{{M_{YZ}}^2\left(x\right)}{2E_{\mathrm{s}}{I}_Z}dx+{\int }_0^l\frac{{F_{Y\mathrm{N}}}^2\left(x\right)}{2E_{\mathrm{s}}A}dx=\frac{{\sigma }_Y^2{\eta }^2{l}^7\mathrm{c}\mathrm{o}{\mathrm{s}}^4\theta }{24E_{\mathrm{s}}{I}_Z}+\frac{{\sigma }_Y^2{\eta }^2{l}^5\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta }{2E_{\mathrm{s}}A}.$

The displacement of the cell along the $Y$ direction, ${\delta }_Y$, is obtained by Castigliano’s Theorem:

11

${\delta }_Y=\frac{\partial U_Y}{\partial F_Y}=\frac{{\sigma }_Y\eta l^5\mathrm{c}\mathrm{o}{\mathrm{s}}^3\theta }{12E_{\mathrm{s}}{I}_Z}+\frac{{\sigma }_Y\eta l^3\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta }{E_{\mathrm{s}}A}.$

The strain along the $Y$ direction, ${\epsilon }_Y$, is as follows:

12

${\epsilon }_Y=\frac{{\delta }_Y}{\alpha l}=\frac{{\sigma }_Y\eta l^2\mathrm{c}\mathrm{o}\mathrm{s}\theta (\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta +{\beta }^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta )}{E_{\mathrm{s}}{\beta }^3\eta l\alpha }.$

The equivalent elastic modulus along the $Y$ direction, $E_Y^{*}$, is expressed according to Hooke’s law:

13

$E_Y^{\mathrm{*}}=E_{\mathrm{s}}\frac{\alpha {\beta }^3}{\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta +\beta \right)\left(\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta +{\beta }^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \right)},\left(\theta <0°\right).$

For hexagonal and quadrilateral cells, the equivalent elastic modulus along the $Y$ direction is written as:

14

$E_Y^{\mathrm{*}}=E_{\mathrm{s}}\frac{{\beta }^3\left(\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\theta \right)}{\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta +\beta \right)\left(\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta +{\beta }^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \right)},\left(\theta \ge 0°\right).$

2.3. Equivalent shear modulus in the $\mathbf{X}\mathbf{Y}$ plane

The displacement of the cell subjected to the shear stress $\tau$ in Fig. 5(a) is obtained by calculating that of point B. Similarly, the vertical wall is regarded as a fixed-end constraint, and the cell is simplified to a symmetric beam with antisymmetric loads, as shown in Fig. 5(b). The load on point B, $F_{XY}$, is expressed by:

15

$F_{XY}=\tau \eta l\left(\frac{1}{2}\alpha l+l\mathrm{s}\mathrm{i}\mathrm{n}\theta \right).$

The displacement caused by the shear deformation, ${\delta }_{XY}$, is obtained by Castigliano’s Theorem:

16

${\delta }_{XY}=\frac{8\mathrm{c}\mathrm{o}{\mathrm{s}}^2\theta ({\alpha }^3+3{\alpha }^2+6\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta +7\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta )F_{XY}{l}^3}{3({\alpha }^3+12{\alpha }^2+24\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta +16\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta )E_{\mathrm{s}}{I}_Z}.$

Fig. 5Concave hexagonal cell subjected to τ and its equivalent model

a) Concave hexagonal cell subjected to $\tau$

b) Equivalent symmetric beam model

The in-plane shear strain ${\gamma }_{XY}$ is as follows:

17

${\gamma }_{XY}=\frac{{\delta }_{XY}}{2l\mathrm{c}\mathrm{o}\mathrm{s}\theta }.$

The equivalent shear modulus $G_{XY}^{*}$ is expressed according to Hooke’s law:

18

$G_{XY}^{\mathrm{*}}=E_{\mathrm{s}}\frac{{\beta }^3\left({\alpha }^3+12{\alpha }^2+24\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta +16\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \right)}{8\alpha \mathrm{c}\mathrm{o}\mathrm{s}\theta \left({\alpha }^3+3{\alpha }^2+6\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta +7\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \right)},\left(\theta <0°\right).$

For hexagonal and quadrilateral cells, the equivalent shear modulus is written as:

19

$G_{XY}^{\mathrm{*}}=E_{\mathrm{s}}\frac{{\beta }^3\left({\alpha }^3+12{\alpha }^2+24\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta +16\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \right)}{8\mathrm{c}\mathrm{o}\mathrm{s}\theta \left({\alpha }^3+3{\alpha }^2+6\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta +7\mathrm{s}\mathrm{i}{\mathrm{n}}^2\theta \right)\left(\alpha +2\mathrm{s}\mathrm{i}\mathrm{n}\theta \right)},\left(\theta \ge 0°\right).$

3. Optimization method of cell structure

The vertical wall length, wall thickness and inner angle of the basic cell are chosen as the design variables:

where $H$ is the vector of design variables.

The elastic modulus is an important index to measure the strength and stiffness of the honeycomb core. Hence, the design objective is to maximize the equivalent elastic and shear moduli:

where ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$ are the weight coefficients, and ${\lambda }_1+{\lambda }_2+{\lambda }_3=1$.

Too thin or too thick cell wall increases the processing difficulty and results into the solid cell, respectively, and the cell will disappear if the inner angle is too great. Besides, the relative density ${\rho }_s^{*}$ is used to characterize the mass in design space. Hence, $\alpha$, $\beta$, $\theta$ and ${\rho }_s^{*}$ are bounded by:

where ${\rho }_{ori}^{*}$ is the relative density of the original cell.

The flow chart of the cell structure optimization is shown in Fig. 6. The multi-objective GA starts the search from a string set of solutions, generating a large number of non-inferior solutions in one optimization thus avoiding the disadvantage that the integrated algorithm is easy to fall into a local optimum. For GA parameters, the crossover probability is 0.7, the mutation probability is 0.02, and the maximum iteration number is 100. To ensure the calculation accuracy, the design variables and optimized results keep 4 decimal places using binary encoding:

where $H_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $H_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are the upper and lower limits of design variables, respectively; $n$ is accurate to the $n$^th decimal place; $t$ is the length of the encoded gene, and is 33.

Fig. 6Flow chart of the cell structure optimization

The optimizer performance is illustrated in Fig. 6(a), and the good convergence is achieved after 19 iterations. The structural parameters are shown in Table 3, and Fig. 7 gives the original and optimized cells.

Fig. 7Optimizer performance and cells

a) Optimizer performance

b) Original and optimized cells

4. Structure performance analysis of original and optimized honeycomb cores

The original and optimized cells are respectively tiled to form the honeycomb cores, and then the finite element models are established with APDL language in ANSYS environment and are meshed by shell93 element, as shown in Fig. 8. The statics and modal analysis are performed with an aim to study the effects of the cell optimization on the displacement, stress, strain, nature frequency and mode shape.

Fig. 8Finite element models of original and optimized honeycomb cores.

a) Before the cell optimization

b) After the cell optimization

Fig. 9Displacement, stress and strain distributions of original and optimized honeycomb cores

a) Displacement before the optimization

b) Displacement after the optimization

c) Stress before the optimization

d) Stress after the optimization

e) Strain before the optimization

f) Strain after the optimization

4.2. Modal analysis

The modal analysis is carried out by Lanczos algorithm, and the frequencies and mode shapes of original and optimized honeycomb cores are obtained. Table 4, Fig. 10 and Fig. 11 give only the first six frequencies and mode shapes due to the greater influence of the lower frequencies on the stability. The first six mode shapes of optimized honeycomb are the same as those of original honeycomb. The second-, third-, fourth- and sixth- order frequencies increase by 2.95 %, 13.38 %, 24.29 % and 2.95 % after the optimization, while the first- and fifth- order frequencies decrease by 2.28 % and 8.33 %.

Table 4Frequencies and mode shapes of original and optimized honeycomb cores

Frequency order	Before the optimization		After the optimization
	Frequency / (Hz)	Mode shape	Frequency / (Hz)	Mode shape
1st	167.39	First-order shimmy vibration	163.57	First-order shimmy vibration
2nd	357.19	First-order torsional vibration	367.73	First-order torsional vibration
3rd	832.18	First-order bending vibration	943.55	First-order bending vibration
4th	1109.41	Second-order bending vibration	1378.84	Second-order bending vibration
5th	1252.38	Second-order shimmy vibration	1148.04	Second-order shimmy vibration
6th	1339.38	Second-order torsional vibration	1378.84	Second-order torsional vibration

Fig. 10Mode shapes of the original honeycomb core

c) First-order

b) Second-order

c) Third-order

d) Fourth-order

e) Fifth-order

f) Sixth-order

Fig. 11Mode shapes of the optimized honeycomb core

c) First-order

b) Second-order

c) Third-order

d) Fourth-order

e) Fifth-order

f) Sixth-order

5. Conclusions

Based on the cantilever and symmetric beam models, the equivalent elastic and shear moduli of the basic cell are derived by Energy Method, respectively. Compared with the method of Ref. [15], the accuracies of the formulas are higher since the influences of tensile and bending deformations of the cell wall are considered. The cell structure parameters are optimized with the consideration of in-plane tensile/compression and shear mechanical properties, and the displacement, stress, strain, frequency and vibration mode of original and optimized honeycomb cores are analyzed. The maximum displacement, stress and strain of the honeycomb core decrease by 50.70 %, 51.69 % and 47.45 % after the cell optimization, thus weakening the stress concentration and deformation and increasing the strength and stiffness. The natural frequencies of optimized honeycomb core mostly increase and are all beyond the excitation frequency range, indicating that the optimized structure is less prone to the resonance.