# Multilayered Equivalent Finite Element Method for Embedded Honeycomb Plates

Authors: Yang Xu; Haohui Wang; Xiaowei Sheng
Journal: Shock and Vibration
 (2018)

---

## Abstract

To investigate the mechanical properties of embedded honeycomb plates with high efficiency and accuracy, a new multilayered equivalent finite element method (FEM) model is proposed. A series of FEM numerical studies (modal analysis, static analysis, and shock spectrum analysis) are performed. The goal is to compare the errors produced by the multilayered equivalent method and by existing equivalent approaches. The obtained results indicate that the proposed model shows good agreement with the original plate. Moreover, based on the new model, a parametric study correlating the microstructure parameters (embedded depth/cell size) to modal frequency is proposed, and a multiparameter equation for frequency and embedded depth/cell size is established to serve as a basis for structural optimization design.


---

## Body

## 1. Introduction

An embedded honeycomb plate is a special honeycomb sandwich structure that is widely used in aerospace engineering owing to its low weight, high specific bending stiffness, strength under distributed loads, and good energy-absorbing capacity.Finite element analysis (FEA) of a honeycomb sandwich panel has been conducted to obtain accurate results for complex problems. However, it requires a long computational time, and the required time increases with problem complexity. Another important problem facing finite element method (FEM) modeling engineers is how to create better designs and to reduce manufacturing costs by predicting the structural performance accurately.To date, many equivalent methods for a honeycomb sandwich plate have been studied. Boudjemai et al. [1] presented a feasible and practical whole plate equivalent method. This approach can only describe macrolevel characteristics and not microlevel ones. Xie et al. [2] investigated the damage and high-speed impact deformation mechanisms in honeycomb sandwich panels. Jeong et al. [3] presented a simple mathematical model and numerical simulations of the hexagonal pattern formation of a honeycomb using the immersed boundary method. Allen [4] developed a honeycomb plate model that ignores in-plane stiffness and bending stiffness. Gibson [5, 6] provided a theoretical equivalent parameter formula for a honeycomb material. Gibson’s formula is widely used owing to its simplicity. Xia and Jin [7] proposed three equivalent methods. These three methods have been used to calculate the natural frequencies of a honeycomb sandwich plate. The computational results show that the three equivalent methods are reliable and practical.This study proposes an analytical method with a new multilayered equivalent FEM model. The feasibility and accuracy of the proposed method were verified by comparing the calculation results obtained using the proposed FEM model with the current FEM results. The honeycomb panel geometry effect was also studied. This study can serve as a guide for optimizing the honeycomb sandwich structure.


## 2. Modeling of Embedded Honeycomb Plates

### 2.1. Embedded Honeycomb Structure

An embedded honeycomb structure consists of two thin face sheets attached to both sides of a lightweight core and embedded parts, as shown in Figure1. This design allows the outer face sheets to carry axial loads, bending moments, and in-plane shears, whereas the core carries normal flexural shears.Figure 1
Embedded honeycomb structure.Honeycomb structures are susceptible to failures due to large normal local stress concentrations because of the nonuniformity of the core/face sheet assembly; therefore, the embedded part is necessary. To maintain homogeneity, the core, skin, and embedded part are typically fabricated using aluminum.


### 2.2. Multilayered Equivalent of Embedded Honeycomb Structure

The sandwich structure [8, 9] can be modeled using the whole entity modeling (Figure 2) and equivalent plate methods (Figure 3).(1)
Whole entity modeling: the skin, core, and embedded part are established using 3D solid elements in FEM software. Though this model has high consistency with the original plate, it requires too much computational time in engineering application.(2)
Equivalent plate: the original plate is equivalent to an isotropic/anisotropic plate. One idealizes the real constitution of the material by considering it continuous. This method can be used to reduce the time spent for the analysis of the honeycomb sandwich structure, and it is greatly advantageous in that it reduces the preprocessor time and computational time. However, the error (modal analysis error, static analysis error, and shock spectrum response analysis error) in the accuracy will become obvious with an increase in the model size.Figure 2
Whole entity modeling.Figure 3
Equivalent plate.These methods, respectively, lack efficiency and accuracy for more complex problems.This study proposes a new multilayered equivalent FEM model to overcome these disadvantages. The principle is that the sandwich structure is equivalent to a laminated material (Figure4). Therefore, it can output the response between layers; this response can be used for analyzing the state of the core [10].Figure 4
Multilayered equivalent.The multilayered equivalent parameters of the plate had to be obtained before performing an analysis using HyperWorks.A theoretical model of a honeycomb cell is established to calculate the equivalent parameters (Figure5). The properties of the honeycomb material under transverse loading are considered, and the homogeneous material is assumed to be uniform in single tensile state (Figure 6); this gives [5–7](1)M=12PLsinθ,(2)P=σ1H+Lsinθb,where b is the panel width, and the deflection of AB based on bending theory is given by(3)ε1=PL3sinθ12ESl=PL3sinθESbt3.Figure 5
Honeycomb cell.Figure 6
Horizontal deformation of honeycomb cell.I
=
1
/
12
b
t
3 is the moment of inertia, and the extension of AB caused by loading P is obtained by(4)ε2=σ1ESL=LPcosθESbt.The equivalent effect becomes(5)δ1=ε1sinθ+ε2cosθLcosθ=Psin2θL3ESbLt3cosθ1+cotθt2L2,(6)δ2=ε2sinθ−ε1cosθH+Lsinθ=−PsinθcosθL3ESbH+Lsinθt31−t2L2.The equivalent Poisson’s ratioμ1 is expressed as follows:(7)μ1=−δ2δ1=cos2θH/L+sinθsinθ·1−t2/L21+cot2θt2/L2.From (2) and (5), the transverse Young’s modulus is expressed as(8)E1=ESt3cosθL3H/L+sinθsin2θ·11+cot2θt2/L2.For analyzing the mechanical properties under longitudinal loading (Figure7), it is also assumed that the equivalent homogeneous material is uniform in single tensile state, unlike transverse loading, and the stretching deformation of BC is considered; therefore,(9)M=12WLsinθ,W=Lσ2sinθ.Figure 7
Longitudinal deformation of honeycomb cell.The deflection of AB based on bending theory is given by(10)ε1=WL3cosθESbt3.The extension of AB caused by loadingW is given by(11)ε2=WLsinθESbt.The elongation of BC is given by(12)ε3=WHESbt.The equivalent strain is then expressed as(13)δ1=−ε1sinθ−ε2cosθLcosθ=−WL3sinθESbt1−t2L2,δ2=ε2sinθ+ε1cosθ+ε3H+Lsinθ=WL4cos2θESbt3H+Lsinθ·1+Hsec2θL+cot2θt2L2.The equivalent Poisson’s ratioμ2 is expressed as(14)μ2=−δ1δ2=H/L+sinθsinθcos2θ·11+Hsec2θ/L+cot2θt2/L2.The longitudinal Young’s modulus is expressed as(15)E2=ESt3H/L+sinθL3cosθ·11+Hsec2θ/L+cot2θt2/L2.Therefore, the material parameters of the honeycomb are as follows:(16)E1=ESt3cosθL3H/L+sinθsin2θ·11+cot2θt2/L2,μ1=−cos2θH/L+sinθsinθ·1−t2/L21+cot2θt2/L2,E2=ESt3H/L+sinθL3cosθ·11+Hsec2θ/L+cot2θt2/L2,μ2=H/L+sinθsinθcos2θ·11+Hsec2θ/L+cot2θt2/L2.For a quantitya that is much lower than 1, the following approximate formula may be used:(17)11+a=1−a,t2/L2 is very small. For the positive hexagon =L, θ=π/6, Equation (16) can be approximately simplified as(18)E1=43ES1−3t2L2t3L3,μ1=1−4t2L2,E2=43ES1−53t2L2t3L3,μ2=1−83t2L2.It is important to consider the stretching deformation of the honeycomb cell (BC) in the aerospace industry because the first few modes of such equipment are usually telescopic and torsional.When the equivalent parameters are input into the multilayered equivalent model, the calculation scale reduces significantly compared with that in whole entity modeling. Compared with the equivalent plate method that also adopts shell element modeling, the proposed model can make the response output among layers greatly improve the accuracy.


## 2.1. Embedded Honeycomb Structure

An embedded honeycomb structure consists of two thin face sheets attached to both sides of a lightweight core and embedded parts, as shown in Figure1. This design allows the outer face sheets to carry axial loads, bending moments, and in-plane shears, whereas the core carries normal flexural shears.Figure 1
Embedded honeycomb structure.Honeycomb structures are susceptible to failures due to large normal local stress concentrations because of the nonuniformity of the core/face sheet assembly; therefore, the embedded part is necessary. To maintain homogeneity, the core, skin, and embedded part are typically fabricated using aluminum.


## 2.2. Multilayered Equivalent of Embedded Honeycomb Structure

The sandwich structure [8, 9] can be modeled using the whole entity modeling (Figure 2) and equivalent plate methods (Figure 3).(1)
Whole entity modeling: the skin, core, and embedded part are established using 3D solid elements in FEM software. Though this model has high consistency with the original plate, it requires too much computational time in engineering application.(2)
Equivalent plate: the original plate is equivalent to an isotropic/anisotropic plate. One idealizes the real constitution of the material by considering it continuous. This method can be used to reduce the time spent for the analysis of the honeycomb sandwich structure, and it is greatly advantageous in that it reduces the preprocessor time and computational time. However, the error (modal analysis error, static analysis error, and shock spectrum response analysis error) in the accuracy will become obvious with an increase in the model size.Figure 2
Whole entity modeling.Figure 3
Equivalent plate.These methods, respectively, lack efficiency and accuracy for more complex problems.This study proposes a new multilayered equivalent FEM model to overcome these disadvantages. The principle is that the sandwich structure is equivalent to a laminated material (Figure4). Therefore, it can output the response between layers; this response can be used for analyzing the state of the core [10].Figure 4
Multilayered equivalent.The multilayered equivalent parameters of the plate had to be obtained before performing an analysis using HyperWorks.A theoretical model of a honeycomb cell is established to calculate the equivalent parameters (Figure5). The properties of the honeycomb material under transverse loading are considered, and the homogeneous material is assumed to be uniform in single tensile state (Figure 6); this gives [5–7](1)M=12PLsinθ,(2)P=σ1H+Lsinθb,where b is the panel width, and the deflection of AB based on bending theory is given by(3)ε1=PL3sinθ12ESl=PL3sinθESbt3.Figure 5
Honeycomb cell.Figure 6
Horizontal deformation of honeycomb cell.I
=
1
/
12
b
t
3 is the moment of inertia, and the extension of AB caused by loading P is obtained by(4)ε2=σ1ESL=LPcosθESbt.The equivalent effect becomes(5)δ1=ε1sinθ+ε2cosθLcosθ=Psin2θL3ESbLt3cosθ1+cotθt2L2,(6)δ2=ε2sinθ−ε1cosθH+Lsinθ=−PsinθcosθL3ESbH+Lsinθt31−t2L2.The equivalent Poisson’s ratioμ1 is expressed as follows:(7)μ1=−δ2δ1=cos2θH/L+sinθsinθ·1−t2/L21+cot2θt2/L2.From (2) and (5), the transverse Young’s modulus is expressed as(8)E1=ESt3cosθL3H/L+sinθsin2θ·11+cot2θt2/L2.For analyzing the mechanical properties under longitudinal loading (Figure7), it is also assumed that the equivalent homogeneous material is uniform in single tensile state, unlike transverse loading, and the stretching deformation of BC is considered; therefore,(9)M=12WLsinθ,W=Lσ2sinθ.Figure 7
Longitudinal deformation of honeycomb cell.The deflection of AB based on bending theory is given by(10)ε1=WL3cosθESbt3.The extension of AB caused by loadingW is given by(11)ε2=WLsinθESbt.The elongation of BC is given by(12)ε3=WHESbt.The equivalent strain is then expressed as(13)δ1=−ε1sinθ−ε2cosθLcosθ=−WL3sinθESbt1−t2L2,δ2=ε2sinθ+ε1cosθ+ε3H+Lsinθ=WL4cos2θESbt3H+Lsinθ·1+Hsec2θL+cot2θt2L2.The equivalent Poisson’s ratioμ2 is expressed as(14)μ2=−δ1δ2=H/L+sinθsinθcos2θ·11+Hsec2θ/L+cot2θt2/L2.The longitudinal Young’s modulus is expressed as(15)E2=ESt3H/L+sinθL3cosθ·11+Hsec2θ/L+cot2θt2/L2.Therefore, the material parameters of the honeycomb are as follows:(16)E1=ESt3cosθL3H/L+sinθsin2θ·11+cot2θt2/L2,μ1=−cos2θH/L+sinθsinθ·1−t2/L21+cot2θt2/L2,E2=ESt3H/L+sinθL3cosθ·11+Hsec2θ/L+cot2θt2/L2,μ2=H/L+sinθsinθcos2θ·11+Hsec2θ/L+cot2θt2/L2.For a quantitya that is much lower than 1, the following approximate formula may be used:(17)11+a=1−a,t2/L2 is very small. For the positive hexagon =L, θ=π/6, Equation (16) can be approximately simplified as(18)E1=43ES1−3t2L2t3L3,μ1=1−4t2L2,E2=43ES1−53t2L2t3L3,μ2=1−83t2L2.It is important to consider the stretching deformation of the honeycomb cell (BC) in the aerospace industry because the first few modes of such equipment are usually telescopic and torsional.When the equivalent parameters are input into the multilayered equivalent model, the calculation scale reduces significantly compared with that in whole entity modeling. Compared with the equivalent plate method that also adopts shell element modeling, the proposed model can make the response output among layers greatly improve the accuracy.


## 3. Examples and Discussion

This section presents the modal analysis, static analysis, and shock spectrum analysis for the embedded honeycomb using three methods: whole entity modeling, equivalent plate, and multilayered equivalent. The whole entity modeling results are used as a benchmark for comparison [11–17].Generally, the epoxy resin adhesive between the skin and the core is important because delaminations (skin-to-core bonding failure) may occur. The adhesive was modeled as a “dummy” layer that is placed directly in contact with the skin and core/part. Embedded honeycomb plates can show several different types of but interacting failure: failure of layers belonging to skins, failure of layers belonging to core/part, and interlaminar failure of skin-to-core bonding. The interlaminar stresses responsible for interlaminar failure cannot be evaluated by employing the classic theory of laminated plates because the latter does not allow the shear stresses at the interface of each layer to be estimated. Moreover, [18, 19] show that the adhesive has no noticeable effects on the final results of tests in a simple loaded state.In this analysis, the adhesive is not considered because its effect is negligible.The plate has length, width, radius, thickness, skin thickness, cell thickness, and cell size of 250 mm, 200 mm, 40 mm, 15 mm, 0.3 mm, 0.04 mm, and 4 mm, respectively, and it is made of aluminum.Constraints:In modal analysis, AB and A1B1 are fixed.In static analysis, A, B, C, D, A1, B1, C1, and D1 are fixed.In shock spectrum analysis, A, B, C, D, A1, B1, C1, and D1 are loading points (Figure 8).Figure 8
Constraints.
(a)
(b)


### 3.1. Modal Analysis

Table1 shows the first five natural frequencies obtained with the three methods used for the modal analysis of the embedded honeycomb plate, and Figure 9 shows the first five mode shapes.Table 1
Natural frequency.
Frequency (Hz)
Method
Whole entity modeling
Equivalent plate
Error
Multilayered equivalent
Error
f
1
95.30
87.09
−8.21
96.77
1.47
f
2
260.72
226.83
−33.89
277.58
16.86
f
3
510.19
432.55
−77.64
522.85
12.66
f
4
570.19
559.41
−10.78
560.16
−10.03
f
5
761.90
707.54
−54.36
811.83
49.93
Maximum % error
—
—
−14.21
—
6.5Figure 9
Various mode shapes. (a) Whole entity modeling, (b) equivalent plate, and (c) multilayered equivalent.
(a)
(b)
(c)According to Table1, the proposed multilayered equivalent method usually produces higher frequencies, indicating a stiffer model. By contrast, the equivalent plate model is consistently more flexible and shows smaller frequencies. For the multilayered equivalent method, it is assumed that the core can resist the transverse shearing deformation and that it has some in-plane stiffness, whereas the top and bottom surface layers cannot resist the shearing deformation but satisfy the Kirchhoff hypothesis. Therefore, the multilayered equivalent is stiffer. The equivalent plate is analyzed using the REDDY low-order shear deformation theory; the honeycomb sandwich panel is replaced with an orthotropic plate of the same stiffness and size in the software. In fact, the equivalent plate is not a completely orthotropic plate; it is somewhere between an isotropic plate and an orthotropic plate. Therefore, the equivalent plate is more flexible [20, 21]. Specifically, the assumption regarding skin-core idealization can involve an incorrect assessment of the stiffness and strength properties of the material itself, and the mechanical properties of the skin-core sandwiches strongly depend on the honeycomb production process.Figure9 shows good agreement between the first five mode shapes. Table 1 shows a comparison between the first five natural frequencies. The error of f1, f2, f3, and f5 of the multilayered equivalent is obviously smaller than that of the equivalent plate. However, the error of f4 is very similar. This could be attributable to the deformation of Mode 4 mainly occurring at the core far away from the embedded part, thereby having little effect on the plate.In addition, errors are more significant at high frequencies. Figure10 shows a good comparison between the results of the multilayered equivalent and the equivalent plate for a large frequency range (0–5000 Hz). The multilayered equivalent clearly well represents the embedded honeycomb plate for modal analysis, and it can be used for complex and large honeycomb structures; this reduces the computing time of the analysis and contains less error.Figure 10
Comparison of models for large frequency range.


### 3.2. Static Analysis

The center node of the plate is selected as the displacement reference point (load isFy = 7.5 g, Fz = 1.2 g). Table 2 shows the results.Table 2
Displacement results.
Parameter (mm)
Method
Whole entity modeling
Equivalent plate
Multilayered equivalent
Y-direction
0.03389
0.03513
0.03315
Z-direction
0.00048
0.00051
0.0004651
Mag direction
0.03391
0.03514
0.03316
Maximum error%
—
6.25
−3.1Table2 shows that the maximum error of the multilayered equivalent (−3.1%) is significantly improved compared with that of the equivalent plate (6.25%).The center node of the plate is selected as the stress reference point; Table3 shows the results.Table 3
Stress result.
Parameter (MPa)
Method
Whole entity modeling
Equivalent plate
Multilayered equivalent
Stress result
1.7
1.54
1.61
Relative error%
—
−9.41
−5.3Table3 shows that the result of the multilayered equivalent is significantly improved compared with that of the equivalent plate; however, the absolute error is still 5.3% because of the difference in boundary conditions in the equivalent process.


### 3.3. Shock Spectrum Response Analysis

The response spectrum of the impulse acceleration is 100–1500 Hz for 80 g and 1500–5000 Hz for 2000 g, and the loading direction is they-direction. Figures 11 and 12, respectively, show the maximum displacement and stress obtained through the shock spectrum analysis.Figure 11
Displacement obtained through shock spectrum analysis. (a) Whole entity modeling, (b) equivalent plate, and (c) multilayered equivalent.
(a)
(b)
(c)Figure 12
Stress obtained through shock spectrum analysis. (a) Whole entity modeling, (b) equivalent plate, and (c) multilayered equivalent.
(a)
(b)
(c)These figures show that the deformation of the plate is slightly different. The displacement deviation near the connections of the embedded part and panel increased because the connections mainly included the effect of the embedded part on the panel, displacement consistency is better, and stress distribution is basically similar.


## 3.1. Modal Analysis

Table1 shows the first five natural frequencies obtained with the three methods used for the modal analysis of the embedded honeycomb plate, and Figure 9 shows the first five mode shapes.Table 1
Natural frequency.
Frequency (Hz)
Method
Whole entity modeling
Equivalent plate
Error
Multilayered equivalent
Error
f
1
95.30
87.09
−8.21
96.77
1.47
f
2
260.72
226.83
−33.89
277.58
16.86
f
3
510.19
432.55
−77.64
522.85
12.66
f
4
570.19
559.41
−10.78
560.16
−10.03
f
5
761.90
707.54
−54.36
811.83
49.93
Maximum % error
—
—
−14.21
—
6.5Figure 9
Various mode shapes. (a) Whole entity modeling, (b) equivalent plate, and (c) multilayered equivalent.
(a)
(b)
(c)According to Table1, the proposed multilayered equivalent method usually produces higher frequencies, indicating a stiffer model. By contrast, the equivalent plate model is consistently more flexible and shows smaller frequencies. For the multilayered equivalent method, it is assumed that the core can resist the transverse shearing deformation and that it has some in-plane stiffness, whereas the top and bottom surface layers cannot resist the shearing deformation but satisfy the Kirchhoff hypothesis. Therefore, the multilayered equivalent is stiffer. The equivalent plate is analyzed using the REDDY low-order shear deformation theory; the honeycomb sandwich panel is replaced with an orthotropic plate of the same stiffness and size in the software. In fact, the equivalent plate is not a completely orthotropic plate; it is somewhere between an isotropic plate and an orthotropic plate. Therefore, the equivalent plate is more flexible [20, 21]. Specifically, the assumption regarding skin-core idealization can involve an incorrect assessment of the stiffness and strength properties of the material itself, and the mechanical properties of the skin-core sandwiches strongly depend on the honeycomb production process.Figure9 shows good agreement between the first five mode shapes. Table 1 shows a comparison between the first five natural frequencies. The error of f1, f2, f3, and f5 of the multilayered equivalent is obviously smaller than that of the equivalent plate. However, the error of f4 is very similar. This could be attributable to the deformation of Mode 4 mainly occurring at the core far away from the embedded part, thereby having little effect on the plate.In addition, errors are more significant at high frequencies. Figure10 shows a good comparison between the results of the multilayered equivalent and the equivalent plate for a large frequency range (0–5000 Hz). The multilayered equivalent clearly well represents the embedded honeycomb plate for modal analysis, and it can be used for complex and large honeycomb structures; this reduces the computing time of the analysis and contains less error.Figure 10
Comparison of models for large frequency range.


## 3.2. Static Analysis

The center node of the plate is selected as the displacement reference point (load isFy = 7.5 g, Fz = 1.2 g). Table 2 shows the results.Table 2
Displacement results.
Parameter (mm)
Method
Whole entity modeling
Equivalent plate
Multilayered equivalent
Y-direction
0.03389
0.03513
0.03315
Z-direction
0.00048
0.00051
0.0004651
Mag direction
0.03391
0.03514
0.03316
Maximum error%
—
6.25
−3.1Table2 shows that the maximum error of the multilayered equivalent (−3.1%) is significantly improved compared with that of the equivalent plate (6.25%).The center node of the plate is selected as the stress reference point; Table3 shows the results.Table 3
Stress result.
Parameter (MPa)
Method
Whole entity modeling
Equivalent plate
Multilayered equivalent
Stress result
1.7
1.54
1.61
Relative error%
—
−9.41
−5.3Table3 shows that the result of the multilayered equivalent is significantly improved compared with that of the equivalent plate; however, the absolute error is still 5.3% because of the difference in boundary conditions in the equivalent process.


## 3.3. Shock Spectrum Response Analysis

The response spectrum of the impulse acceleration is 100–1500 Hz for 80 g and 1500–5000 Hz for 2000 g, and the loading direction is they-direction. Figures 11 and 12, respectively, show the maximum displacement and stress obtained through the shock spectrum analysis.Figure 11
Displacement obtained through shock spectrum analysis. (a) Whole entity modeling, (b) equivalent plate, and (c) multilayered equivalent.
(a)
(b)
(c)Figure 12
Stress obtained through shock spectrum analysis. (a) Whole entity modeling, (b) equivalent plate, and (c) multilayered equivalent.
(a)
(b)
(c)These figures show that the deformation of the plate is slightly different. The displacement deviation near the connections of the embedded part and panel increased because the connections mainly included the effect of the embedded part on the panel, displacement consistency is better, and stress distribution is basically similar.


## 4. Embedded Honeycomb Structure Geometry Effect

Many design parameters can affect the natural frequency of the honeycomb plate, including the skin thickness, core thickness, cell size, and embedded depth.Most studies on the optimization of honeycomb panels have rarely used the cell size or embedded depth as parameters. Both are geometry parameters, and changing them involves redesigning the honeycomb core; this is equivalent to remodeling every iteration, which existing FEA software can hardly do [22–25].Based on the multilayered equivalent FEM model, equivalent material parameters are used as variables. The equivalent parameter formula is linear and reversible [26], and each change in cell size/embedded depth leads to a change in the material parameters; therefore, the relationship between cell size/embedded depth and natural frequency can be understood indirectly.The following subsection discusses the effect of cell size and embedded depth on the natural frequency of an embedded honeycomb plate.


### 4.1. Effect of Cell Size

This section discusses the influence of cell size on natural frequency. The frequency and multiparameter equations can be expressed as(19)fi=fh,l,t,tfi=1,2,3,4,5,where h is the core height; l, the cell size; t, the embedded depth; and tf, the skin thickness.The following equations are obtained:(20)dfi=∂fi∂hdh+∂fi∂ldl+∂fi∂tdt+∂fi∂tfdtf.Equation (20) can be approximately expressed as(21)Δfi=∂fi∂hΔh+∂fi∂lΔl+∂fi∂tΔt+∂fi∂tfΔtf.Equation (21) shows that the change in the natural frequency results from a change in one or more structural parameters.Cell sizes of 1–10.5 mm are tested while keeping other parameters constant. The calculations are imported into Matlab. Figure13 shows the relationship between the first five frequencies and cell size.Figure 13
Relationship between frequency and cell size.The determination of the constant parameter depends on the relation between the parameter and the frequency in a certain region. For data points in the same range, more than one function relationship may exist. However, these relationships represent the same parameter meaning in this region. Therefore, for determining the parameter range derived from the data points to fit the equation, we only need to choose a simple and convenient continuity equation, that is, a polynomial form (Table4).Table 4
Function relationship.
Function
Mode
Constant
a
1
a
2
a
3
a
4
a
5
a
6
f
i
=
a
1
l
5
+
a
2
l
4
+
a
3
l
3
+
a
4
l
2
+
a
5
l
+
a
6
i
=
1
−
4.401
×
10
−
5
0.001478
−
0.01947
0.1348
−
1.123
99.13
i
=
2
−
0.0004448
0.01585
−
0.2314
1.936
−
14.14
284.1
i
=
3
0.003191
−
0.1371
2.21
−
16.14
37.61
523.2
i
=
4
−
0.003897
0.1615
−
2.563
19.38
−
69.4
640.9
i
=
5
−
0.002419
0.08857
−
1.348
11.86
−
82.62
983.8Figure13 shows the frequency variation with different cell sizes of the plate for the first five modes. It is clear that the frequencies decrease with an increase in the cell size for modes 3, 4, and 5 and show a small increase for modes 1 and 2.The cell size does not have a major influence on modes 1 and 2. From mode 3, the frequency change gradually increases. In mode 5, the initial change is the most obvious, and it tends to be steady with an increase in the cell size.The influence of cell size on modes 3, 4, and 5 is greater than that on modes 1 and 2 [27–29]. This is because modes 1 and 2 mainly show bending deformation, and modes 3, 4, and 5 are telescopic and torsional (Figure 9). As noted in Section 2.2, stretching deformation (BC) of cells has greater influence on modes 3, 4, and 5.


### 4.2. Effect of Embedded Depth

This section discusses the influence of embedded depth on natural frequency.By using (19)–(21) and as discussed in Section 4.1, embedded depths of 0.72–14.5 mm are used while keeping other parameters constant. Figure 14 and Table 5 show the results.Figure 14
Relationship between frequency and embedded depth.Table 5
Function relationship.
Function
Mode
Constant
a
1
a
2
a
3
a
4
a
5
a
6
f
i
=
a
1
t
5
+
a
2
t
4
+
a
3
t
3
+
a
4
t
2
+
a
5
t
+
a
6
i
=
1
0.0001794
−0.01007
0.1389
1.461
−41.39
569.4
i
=
2
6.997 × 10−5
−0.002908
−0.02396
2.896
−44.97
593.7
i
=
3
0.0002212
−0.01331
0.2323
0.251
−35.12
640.7
i
=
4
0.0001616
−0.009187
0.1354
0.997
−33.39
650.7
i
=
5
0.0001318
−0.006867
0.06889
1.894
−39.09
678.5Figure14 shows the frequency variation with different embedded depths of the plate for the first five modes. The frequencies clearly decrease (increase) with an increase in the embedded depth for modes 1, 2, and 4. The frequencies of modes 3 and 5 first increase and then decrease, and the inflection point is located at ∼7.5 mm depth, indicating a certain mutation between the embedded depth and the plate stiffness [11, 30–33].A theoretical model of an embedded section with unit width is established (Figure15). Here, M is the upper skin thickness + embedded depth, total plate thickness is 15 mm, and lower skin thickness is 0.3 mm.Figure 15
Theoretical model.This model is equivalent to honeycomb panels with unequal upper and lower skin thicknesses. Here, M∈[0.3, 14.7].M = 0.3 mm indicates a honeycomb panel with equal upper and lower skin thicknesses. M = 14.7 mm indicates a fully penetrating embedded entity. From Figure 9, modes 3 and 5 mainly show bending deformation, following which bending stiffness is important. The bending stiffness of the plate [34, 35] is expressed as(22)D=EI,where D is the bending stiffness of the plate; E, the elastic modulus; and I, the moment of inertia.(23)I=IM+Id,where IM and Id are the inertial moments of the neutral axis, and they are, respectively, defined as(24)IM=Iyc+a2A1,Id=Iyc′+b2A2,where a and b are the distances from yc and yc′ to the y-axis, respectively, and A1 and A2 are cross-sectional areas. These equations are expressed as(25)dIM=H−M22Mdx+M312dx,(26)dId=0.3×147202dx+0.3312dx,(27)dI=H−M22M+M312dx+0.3×147202+0.3312dx.Because the embedded part of the plate is nonuniform, the value varies depending on the selected part. However, Equation (27) shows that although the integral value is different, it is a constant integral to x. Then, constants G1 and G2 are used:(28)I=G1H−M22M+M312+0.3×147202+0.3312.Equations (28) and (22) are combined to give(29)D=EG1H−M22M+M312+G3,(30)D′=EG1M−H22.Equation (27) shows that the rate of change of D within the value range is extreme (M = 7.5 mm); M = 7.5 mm is exactly where the neutral layer of the plate is located. This illustrates that when the mode is dominated by bending deformation, simply increasing the embedded depth does not stabilize the structure. When the embedded depth is less than the neutral layer of the plate, the plate stiffness increases with the embedded depth. With an increase in depth, the plate stiffness tends to decrease as the embedded depth increases when the embedded depth is greater than the neutral layer.


## 4.1. Effect of Cell Size

This section discusses the influence of cell size on natural frequency. The frequency and multiparameter equations can be expressed as(19)fi=fh,l,t,tfi=1,2,3,4,5,where h is the core height; l, the cell size; t, the embedded depth; and tf, the skin thickness.The following equations are obtained:(20)dfi=∂fi∂hdh+∂fi∂ldl+∂fi∂tdt+∂fi∂tfdtf.Equation (20) can be approximately expressed as(21)Δfi=∂fi∂hΔh+∂fi∂lΔl+∂fi∂tΔt+∂fi∂tfΔtf.Equation (21) shows that the change in the natural frequency results from a change in one or more structural parameters.Cell sizes of 1–10.5 mm are tested while keeping other parameters constant. The calculations are imported into Matlab. Figure13 shows the relationship between the first five frequencies and cell size.Figure 13
Relationship between frequency and cell size.The determination of the constant parameter depends on the relation between the parameter and the frequency in a certain region. For data points in the same range, more than one function relationship may exist. However, these relationships represent the same parameter meaning in this region. Therefore, for determining the parameter range derived from the data points to fit the equation, we only need to choose a simple and convenient continuity equation, that is, a polynomial form (Table4).Table 4
Function relationship.
Function
Mode
Constant
a
1
a
2
a
3
a
4
a
5
a
6
f
i
=
a
1
l
5
+
a
2
l
4
+
a
3
l
3
+
a
4
l
2
+
a
5
l
+
a
6
i
=
1
−
4.401
×
10
−
5
0.001478
−
0.01947
0.1348
−
1.123
99.13
i
=
2
−
0.0004448
0.01585
−
0.2314
1.936
−
14.14
284.1
i
=
3
0.003191
−
0.1371
2.21
−
16.14
37.61
523.2
i
=
4
−
0.003897
0.1615
−
2.563
19.38
−
69.4
640.9
i
=
5
−
0.002419
0.08857
−
1.348
11.86
−
82.62
983.8Figure13 shows the frequency variation with different cell sizes of the plate for the first five modes. It is clear that the frequencies decrease with an increase in the cell size for modes 3, 4, and 5 and show a small increase for modes 1 and 2.The cell size does not have a major influence on modes 1 and 2. From mode 3, the frequency change gradually increases. In mode 5, the initial change is the most obvious, and it tends to be steady with an increase in the cell size.The influence of cell size on modes 3, 4, and 5 is greater than that on modes 1 and 2 [27–29]. This is because modes 1 and 2 mainly show bending deformation, and modes 3, 4, and 5 are telescopic and torsional (Figure 9). As noted in Section 2.2, stretching deformation (BC) of cells has greater influence on modes 3, 4, and 5.


## 4.2. Effect of Embedded Depth

This section discusses the influence of embedded depth on natural frequency.By using (19)–(21) and as discussed in Section 4.1, embedded depths of 0.72–14.5 mm are used while keeping other parameters constant. Figure 14 and Table 5 show the results.Figure 14
Relationship between frequency and embedded depth.Table 5
Function relationship.
Function
Mode
Constant
a
1
a
2
a
3
a
4
a
5
a
6
f
i
=
a
1
t
5
+
a
2
t
4
+
a
3
t
3
+
a
4
t
2
+
a
5
t
+
a
6
i
=
1
0.0001794
−0.01007
0.1389
1.461
−41.39
569.4
i
=
2
6.997 × 10−5
−0.002908
−0.02396
2.896
−44.97
593.7
i
=
3
0.0002212
−0.01331
0.2323
0.251
−35.12
640.7
i
=
4
0.0001616
−0.009187
0.1354
0.997
−33.39
650.7
i
=
5
0.0001318
−0.006867
0.06889
1.894
−39.09
678.5Figure14 shows the frequency variation with different embedded depths of the plate for the first five modes. The frequencies clearly decrease (increase) with an increase in the embedded depth for modes 1, 2, and 4. The frequencies of modes 3 and 5 first increase and then decrease, and the inflection point is located at ∼7.5 mm depth, indicating a certain mutation between the embedded depth and the plate stiffness [11, 30–33].A theoretical model of an embedded section with unit width is established (Figure15). Here, M is the upper skin thickness + embedded depth, total plate thickness is 15 mm, and lower skin thickness is 0.3 mm.Figure 15
Theoretical model.This model is equivalent to honeycomb panels with unequal upper and lower skin thicknesses. Here, M∈[0.3, 14.7].M = 0.3 mm indicates a honeycomb panel with equal upper and lower skin thicknesses. M = 14.7 mm indicates a fully penetrating embedded entity. From Figure 9, modes 3 and 5 mainly show bending deformation, following which bending stiffness is important. The bending stiffness of the plate [34, 35] is expressed as(22)D=EI,where D is the bending stiffness of the plate; E, the elastic modulus; and I, the moment of inertia.(23)I=IM+Id,where IM and Id are the inertial moments of the neutral axis, and they are, respectively, defined as(24)IM=Iyc+a2A1,Id=Iyc′+b2A2,where a and b are the distances from yc and yc′ to the y-axis, respectively, and A1 and A2 are cross-sectional areas. These equations are expressed as(25)dIM=H−M22Mdx+M312dx,(26)dId=0.3×147202dx+0.3312dx,(27)dI=H−M22M+M312dx+0.3×147202+0.3312dx.Because the embedded part of the plate is nonuniform, the value varies depending on the selected part. However, Equation (27) shows that although the integral value is different, it is a constant integral to x. Then, constants G1 and G2 are used:(28)I=G1H−M22M+M312+0.3×147202+0.3312.Equations (28) and (22) are combined to give(29)D=EG1H−M22M+M312+G3,(30)D′=EG1M−H22.Equation (27) shows that the rate of change of D within the value range is extreme (M = 7.5 mm); M = 7.5 mm is exactly where the neutral layer of the plate is located. This illustrates that when the mode is dominated by bending deformation, simply increasing the embedded depth does not stabilize the structure. When the embedded depth is less than the neutral layer of the plate, the plate stiffness increases with the embedded depth. With an increase in depth, the plate stiffness tends to decrease as the embedded depth increases when the embedded depth is greater than the neutral layer.


## 5. Conclusion

(1)
Compared with the current modeling methods, the results of modal analysis, static analysis, and shock spectrum analysis obtained by the multilayered equivalent model considering the longitudinal stretching deformation of the honeycomb cell can maximize the accuracy while ensuring efficiency.
(2)
The multilayered equivalent model shows better fitting of the constraints during actual working while ensuring calculation efficiency.
(3)
Based on the multilayered equivalent model, the function relationship between the cell size/embedded depth and frequency is established to intuitively reflect the influence of parameters on the natural frequency. This study can serve as a guide for optimizing the honeycomb sandwich structure.A potential limitation of the proposed method is that the equivalent parameters are identified in a simple and idealized state. This method cannot be used for analyzing problems related to stress concentration and bonding failure (delaminations).



---

*Source: 5914137-2018-09-19.xml*