# Vibroacoustic flexural properties of symmetric honeycomb sandwich panels with composite faces

Laurent Guillaumie

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laurent Guillaumie. Vibroacoustic flexural properties of symmetric honeycomb sandwich panels with composite faces. Journal of Sound and Vibration, 2015, 343, pp.71-103. 10.1016/j.jsv.2014.12.026. hal-01225339

HAL Id: hal-01225339

https://hal.science/hal-01225339

Submitted on 8 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Vibroacoustic flexural properties of symmetric honeycomb sandwich panels with composite faces

Laurent Guillaumie

ONERA – The French Aerospace Lab, BP 72, F-92322 Châtillon, France

The vibroacoustic bending properties of honeycomb sandwich panels with composite faces are studied from the wavenumber modulus to the mechanical impedance, passing through the modal density. Numerical results extracted from finite element software computations are compared with analytical results. In both cases, the homogenization method is used to calculate the global properties of the sandwich panel. Since faces are made of composite material, the classical laminate theory serves as reference. With particular conditions used in the application for symmetric panels, the original orthotropic mechanical properties can be reduced simply to three parameters commonly used in vibroacoustic characterizations. These three parameters are the mass per unit area, the bending rigidity and the out-of-plane shear rigidity. They simultaneously govern the wavenumber modulus, the modal frequencies, the modal density and the mechanical impedance. For all of these vibroacoustic characterizations, a special frequency called the transition frequency separates two domains. In the first domain, below the transition frequency or for low frequencies, the orthotropic sandwich panel has a classical isotropic plate behavior. In the second domain, above the transition frequency or for high frequencies, the out-of-plane shear rigidity is very significant and changes the behavior.

However, the results discussed are only valid up to a certain frequency which is determined by the thickness and out-of-plane shear stiffness of the honeycomb core, the thickness and the bending stiffness of the laminated face sheets and then the mass per unit area and bending stiffness of the total sandwich structure. All these parameters influence the final choice of model and simplifications presented.

Experimental measurements of the bending wavenumber modulus and modal frequencies for our own application were carried out. In the vibroacoustic domain, the critical frequency is also an important frequency. It again depends on the mass per unit area, the bending rigidity and the out-of-plane shear rigidity. The experimental and numerical results of the article are reasonably in agreement with the analytical formula. They all confirm the changes in frequency through different boundary conditions around the panel.

The analytical modal frequencies of rectangular sandwich panels with transverse shear, under simply supported boundary conditions are well known, but under free boundary conditions it is more difficult to predict them. For experiments, however, these latter conditions are the most common. We present, in this paper, an analytical formula that we have developed for the modal frequencies of such a panel under free boundary conditions. All parameters being controlled, it is possible from dynamic measurements and with a special process to identify some honeycomb and composite mechanical properties.
Nomenclature

a, b dimensions of the panel, in the reference x and y directions

A area of the panel

A^M membrane matrix of the panel

A^hom homogenized membrane matrix A^M

A_{a\beta} separate elements of the matrix A^M

c_p phase velocity

c_g group velocity

C^S transverse shear matrix of the panel

C^hom homogenized transverse shear matrix C^S

C_{a\beta} separate elements of the matrix C^S

C_{11}^{hom}, C_{22}^{hom}, C_{12}^{hom} separate homogenized transverse shear elements of the matrix C^{hom}

C_{\phi,ij}^S or C_{\phi,m,n}^S correction factor with transverse shear S of the modal frequencies, under boundary conditions φ

D^F bending matrix of the orthotropic panel

D^hom homogenized bending matrix D^F

D simple common element for bending matrix D^F of the orthotropic panel, with isotropic plate behavior

D_{a\beta} separate elements of the matrix D^F

D_{a\beta}^S elements of the matrix D^F in column

e face sheet thickness of the sandwich plate

e_k thickness of ply k of a multilayer

e_p ply composite thickness

e_{UD} unidirectional composite thickness

E isotropic material Young modulus

E_{L}^{H}, E_{T}^{H}, E_{T}^{UD} honeycomb Young moduli along material axes L, T and U

E_{1}^{UD}, E_{2}^{UD}, E_{3}^{UD} unidirectional composite Young moduli along material axes 1, 2 and 3

f frequency

f_trans transition frequency of the panel

f_cr critical frequency of the panel

f_SSP, i,j modal frequencies of the panel, under simply supported boundary conditions

f_{\phi,ij}^S or f_{\phi,m,n}^S modal frequencies (i, j) or (m, n) of the panel with transverse shear S, under boundary conditions φ

G isotropic material shear modulus

G_{xz}^c, G_{xz}^c transverse shear moduli of the core in reference planes xz and yz

G_g^c geometric transverse shear modulus of the core

G_{12}^{UD} unidirectional composite shear modulus in material plane 12

G_{13}^k, G_{23}^k transverse shear moduli of ply k of a multilayer in material planes 13 and 23

G_{13}^{UD}, G_{23}^{UD} unidirectional composite transverse shear moduli in material planes 13 and 23

G_{LT}^H honeycomb shear modulus in material plane LT

G_{LT'}^H, G_{TT'}^H honeycomb transverse shear moduli in material planes LT' and TT'

G_{xz}^H, G_{yz}^H honeycomb transverse shear moduli in reference planes xz and yz

G_g^H geometric transverse shear modulus of the honeycomb

G_{xz}^k, G_{yz}^k transverse shear moduli of ply k of a multilayer in reference planes xz and yz

G_{xzeq}, G_{yzeq} three-layer transverse shear moduli in reference planes xz and yz

h_H honeycomb thickness

h_S sandwich thickness

i, j half-stationary sine wavenumber in the reference x and y directions

I_1, I_3 first- and third-order mass moment of inertia of the panel

k_c isotropic transverse shear shape coefficient

k_{xz}, k_{yz} three-layer transverse shear shape coefficients in reference planes xz and yz

k_{11}, k_{22}, k_{12} multilayer composite transverse shear shape coefficients in reference planes xz and yz

k(ω) or k(f) bending wavenumber in the panel, at circular frequency ω or at frequency f

k_{ac}(ω) or k_{ac}(f) sound wavenumber, at circular frequency ω or at frequency f

k_x, k_y wavenumbers in the reference directions x and y of the panel

m, n number of nodal lines in the reference x and y directions

m_s mass per unit area of the panel

S reference transverse shear stiffness of the panel

S_{xz}, S_{yz} transverse shear stiffness of the panel in planes xz and yz

n(ω) or n(f) bending modal density of the panel, at circular frequency ω or at frequency f

Q_i reduced stiffness matrix of the isotropic material

Q_{a\beta}^I separate elements of the matrix Q_i

Q_{a\beta}^{I} elements of the matrix Q_i in column

Q_k^*, Q_k material and reference reduced stiffness matrices of ply k of a multilayer

Q_{a\beta}^{*k}, Q_{a\beta}^{k} separate elements of the matrices Q_k^* and Q_k
{Q_{a\beta}^{*}}_k, \{Q_{a\beta}^{*}\}_k elements of the matrices Q_k^* and Q_k in column

Q_H^*, Q_H material and reference reduced stiffness matrices of the honeycomb

Q_{a\beta}^{*H}, Q_{a\beta}^{H} separate elements of the matrices Q_H^* and Q_H

Q_{a\beta}^{*} \{Q_{a\beta}^{*}\}_{H'} elements of the matrices Q_H^* and Q_H in column

Q_{UD}^*, Q_{UD} material and reference reduced stiffness matrices of the unidirectional composite

Q_{a\beta}^{*UD}, Q_{a\beta}^{UD} separate elements of the matrix Q_{UD}^* and Q_{UD}

Q_{a\beta}^{*} \{UD\}, \{Q_{a\beta}^{*}\}_{UD} elements of the matrices Q_{UD}^* and Q_{UD}^* in column

V/F(ω) or V/F(f) mobility of the panel at circular frequency ω or at frequency f

w, \ddot{w} normal displacement and acceleration of the panel
$Z_\infty(\omega)$ or $Z_\infty(f)$ infinite mechanical impedance of the panel, at circular frequency $\omega$ or at frequency $f$

Greek symbols

$\bar{\delta}(f)$ mean frequency interval

$\varphi$ plate boundary conditions (FP: free plate; SSP: simply supported plate)

$\nu$ isotropic material Poisson coefficient

$\nu_{LT}^H, \nu_{TL}^H, \nu_{LT}^H, \nu_{TL}^H, \nu_{TT}^H, \nu_{TT}^H$ honeycomb Poisson coefficients

$\nu_{12}^{UD}, \nu_{21}^{UD}, \nu_{13}^{UD}, \nu_{31}^{UD}, \nu_{23}^{UD}, \nu_{32}^{UD}$ unidirectional composite Poisson coefficients

$\theta_k$ rotation angle of ply $k$ of a multilayer

$\rho$ isotropic material density

$\rho_k$ density of ply $k$ of a multilayer

$\rho_H$ honeycomb density

$\rho_{UD}$ unidirectional composite density

$\omega$ circular frequency ($2\pi f$)

Special functions

$[G_{x_{\varphi,m}}, G_{y_{\varphi,n}}, H_{x_{\varphi,m}}, H_{y_{\varphi,n}}, j_{x_{\varphi,m}}, j_{y_{\varphi,n}}]$ Warburton analytical coefficients for an isotropic plate, as a function of $\varphi$, $m$ and $n$

$[G_{1_{\varphi,j}}, G_{2_{\varphi,j}}, H_{1_{\varphi,j}}, H_{2_{\varphi,j}}, j_{1_{\varphi,j}}, j_{2_{\varphi,j}}]$ Blevins analytical coefficients for an orthotropic plate, as a function of $\varphi$, $i$ and $j$

Cartesian coordinate systems

(1, 2, 3) unidirectional composite axes

$(L, T, T')$ honeycomb axes

$(x, y, z)$ panel reference axes

# 1. Introduction

Sandwich panels are being increasingly used for various industrial applications. They are composed of a relatively soft and thick core made of honeycomb or foam, or both together, lying between two face sheets, which can either be made of an isotropic material or be a stack of laminated composite plies. With their high strength-to-weight ratios, the mechanical and acoustic properties of such honeycomb panels are often difficult to identify. The comparison between numerical model results and dynamical experimental characterizations often shows great discrepancies.

The study of the vibroacoustic bending properties of honeycomb sandwich panels with isotropic face sheets has received much interest in the literature for quite a long time. These panels are generally regarded as a simple uniform three-layer system. These homogenized panels are analyzed within the scope of structural dynamics, rarely in terms of the eigenfrequencies and eigenmodes, as in the work of Wilkinson [1], but more frequently for the statistical energy analysis (SEA Ref. Lyon and Dejong [2]) in terms of wavenumber modulus, modal density and mechanical impedance as in the work of Wilkinson [1], Erickson [3], Mead [4], or Clarkson and Ranky [5]. The differential equation governing the bending displacement uses mainly three parameters: the mass per unit area, the common bending rigidity in the x and y directions and the out-of-plane shear rigidity. When the honeycomb sandwich panels are manufactured with composite skins such as multilayer laminates, the identification of the mechanical properties of the panel with the classical laminate theory (CLPT) quickly becomes complex. To simplify the model, the configuration is rather a symmetric laminate stack so that bending and membrane modes are decoupled. Still for SEA, the authors Renji et al. [6] introduce the orthotropic bending matrix in both directions and include the core transverse shear stiffness in a new analytical model. The wavenumber modulus and modal density are once again studied, detailed and identified, but not the eigenfrequencies. When the eigenfrequencies of orthotropic plates are analyzed, such as by Blevins [7], they are without the transverse shear stiffness term. The sensitivity of the out-of-plane shear relative to the honeycomb thickness and the increasing number of parameters associated with composite laminates must be fully mastered. While there are still large gaps between finite element model results and dynamical experiments, to update the global laminate models with the measured bending mode frequencies, it is necessary to develop the complex method such as the 'Resoanalyzer' or the 'multi-model Resoanalyzer' procedure presented by Lauwagie et al. [8]. However, these methods are only oriented for multilayers of composites without honeycomb, either without transverse shear. Over the last years, Yu and Cleghorn [9] have shown the difficulty of taking into account the out-of-plane shear of symmetric honeycomb panels. To obtain an exact analytical solution, the transverse shear varying with the thickness, we must identify whether the shear deformation is the first- or third-order plate theory (FSDT or TSDT), respectively also called "Reissner-Mindlin" or "Valsov-Levinson-Reddy". The analytical solutions are given for both kinds of shear and for a simply supported rectangular symmetric honeycomb panel, but not under free boundary conditions. Moreover, the chosen models of the plate are tested with skin and honeycomb made of aluminum.

In 1959, Kurtze and Watters [10] published a by now classic paper. The aim of the paper was the development of a simple model for the prediction of sound transmission through sandwich panels. The starting point was a model describing the bending of the panel based on the mechanical impedance of each layer. The laminated face sheets are described as thin skins. Shear effects in the core are considered. The bending stiffness of the sandwich structure is found to vary between two limits; in the low frequency region the bending stiffness is determined by the static properties of the entire structure, the high frequency asymptote is determined by the bending stiffness of the laminated faces. Other models have later confirmed the findings of Kurtze and Watters. During the last decade it has been found that 6th order differential equations derived
based on Hamilton's principle well describe the dynamic and acoustic properties of standard sandwich structures in the entire audible frequency range.

The previous literature showed that the analytical formula for the eigenfrequencies of a symmetric honeycomb sandwich panel with composite face sheets, under simply supported or free boundary conditions, has never been suggested. The main objective of this paper is to propose a formula that takes into account the transverse shear of the plate for two boundary conditions all around the plate (simply supported or free). The formulation that has been developed is however only valid for thin or moderately thick plates. The presentation consists of two parts.

**Part 1:** Firstly, a reminder of the properties of the honeycomb sandwich with isotropic skins is presented (Section 2): wavenumber modulus, modal frequencies, modal density and mechanical impedance. Thereafter we examine the theoretical aspects of a sandwich with composite faces (Section 3). Then, we study a particular application to a symmetric sandwich with composite faces (Section 4) and the technical aspects, in order to compare (Section 5) analytical formulations and numerical computations performed with the finite element code MSC/NASTRAN to experimental results. So far, we make sure that, in the case of our application, all experimental, analytical and numerical SEA results presented herein are consistent with current observations made by the scientific community [5,6]. It is very important to control this set before presenting our proposal.

Here, the wavenumber modulus is not only compared to the experimental and analytical results obtained from a set of rectangular sandwich panels tested in various research projects, but also validated with numerical computations performed using a finite element model and several different boundary conditions applied around the plates. This is not a common presentation. A complementary work provides an analytical formulation of the mechanical impedance, directly as a function of the mass per unit area, the bending rigidity and the transverse shear. In this study, the global analytical transverse shear modulus of the used panels results in an important synthesis of various formulas coming from some authors, such as Wilkinson [1], Allen cited by Erickson [3], Gay [11] and Berthelot [12,13]. The analytical evaluation of the total transverse shear is controlled and is consistent with the results of the numerical model.

**Part 2:** The second part represented by Section 6 explains the original development and choice of the suggested analytical formula for the eigenfrequencies of an orthotropic panel with transverse shear. Inversely from this formula, through experimental eigenfrequencies extracted from a modal analysis and knowing the transverse shear stiffness, we can estimate all elements of the bending matrix with a particular process. The transverse shear stiffness needs to be known for example from the wavenumber modulus experiments, either from the modal density or mechanical impedance measurement. If the faces are made up of a stack of unidirectional composite plies such as [0°,45°, 90°,-45°]s, the in-plane Young and shear moduli can finally be extracted from the elements of the sandwich bending matrix (e.g. $E_{1}^{UD}$, $E_{2}^{UD}$, $G_{12}^{UD}$ and $\nu_{12}^{UD}$, see application: Appendix A). To our knowledge, this is an advance to obtain access to the properties of a unidirectional composite, covering a large frequency range.

The analytical formula for the eigenfrequencies of a rectangular orthotropic panel is presented, using either anterior Blevins coefficients [7], or Warburton coefficients [14] and a particular factor varying with the boundary conditions, called correction factor of transverse shear. All of these extensions are finally applied to our case studied: rectangular symmetric honeycomb sandwich plates with composite faces. To test whether or not the formula yields an accurate prediction, the eigenfrequencies and ranking of eigenmodes are simultaneously observed with the modal analysis results of two finite element models (FEM) and some experiments. It is important to compare analytical and numerical eigenfrequencies with the same classification of the mode ranks. This analysis is rarely presented. This new work is complementary to the developments by Yu and Cleghorn [9] by the extension of their application of a honeycomb sandwich with isotropic face sheets to composite face sheets.

## 2. Symmetric honeycomb sandwich with isotropic faces

In 1967, Wilkinson [1] worked on thin sandwich shells or orthotropic plates with a global transverse shear modulus $G^c$ for the core. Then, Erickson [3] presented modal densities of sandwich panels (curved or flat), in which he used two transverse shear moduli ($G_{xz}^c$, $G_{yz}^c$) for the orthotropic core, in the xz and yz planes. For a multilayer with a core thickness $h_H$ and two equal thickness isotropic face sheets e, the two transverse shear stiffness of the sandwich plate ($S_{xz}$, $S_{yz}$) are based on the Allen formula as

$$S_{xz} = h_H \left(1 + \frac{e}{h_H}\right)^2 G_{xz}^c \quad \text{and} \quad S_{yz} = h_H \left(1 + \frac{e}{h_H}\right)^2 G_{yz}^c \qquad (1)$$

$S_{xz}$ and $S_{yz}$ are used together to define a global (or geometrical) transverse shear stiffness S for the entire plate:

$$S = (S_{xz} \cdot S_{yz})^{1/2} = h_H \left(1 + \frac{e}{h_H}\right)^2 G_g^c \quad \text{with} \quad G_g^c = (G_{xz}^c \cdot G_{yz}^c)^{1/2}. \qquad (2)$$

Later, Mead [4] studied sandwich plates with a honeycomb core and unequal isotropic faces. In this case, all terms of the bending stiffness matrix plate were equal to the same parameter D. The differential equation of motion of such a rectangular three-layer plate is defined taking into account the mass per unit area $m_s$, D and S. The model is normally described by a sixth order differential equation. In 1983, Clarkson and Ranky [5] developed, from this last equation, a simplified differential equation for the displacement with one condition: the face isotropic sheets must be very thin compared to the core
<figure><img src="image_3.png" /><figcaption>Fig. 1. Symmetric honeycomb sandwich with isotropic faces.</figcaption></figure>

thickness. The original model is thus reduced in a 4th order differential equation. This equation can also be used for equal isotropic faces, as in the case of Fig. 1. The honeycomb sandwich is symmetric.

## 2.1. Differential equation of motion

The differential equation for displacement w in the direction z, developed by Clarkson and Ranky is

$$ \underbrace{D\nabla^4 w + m_s \ddot{w}}_{\text{honeycomb sandwich with isotropic faces and with transverse shear}} - m_s \frac{D}{S} \nabla^2 \ddot{w} = 0, \quad (3) $$

with

$$ \nabla^4 = \left( \frac{\partial^4}{\partial x^4} + 2 \frac{\partial^4}{\partial x^2 \partial y^2} + \frac{\partial^4}{\partial y^4} \right), \quad \nabla^2 = \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) \quad \text{and} \quad \ddot{w} = \frac{\partial^2 w}{\partial t^2}. \quad (4) $$

## 2.2. Mass per unit area

The mass per unit area is

$$ m_s = 2\rho e + \rho_H h_H. \quad (5) $$

## 2.3. Bending stiffness

Skins are made of isotropic material. The material reduced stiffness matrix of the face $Q_I$ is defined with $E$, $\nu$ and $G$. The core is an orthotropic honeycomb. The material $L$, $T$ and $T'$ axes of the honeycomb are respectively parallel to the reference $x$, $y$ and $z$ axes of the plate. Assuming that the material and reference axes are the same, the reference reduced stiffness matrix of honeycomb $Q_H$, along axes $x$ and $y$, is considered identical to the material reduced stiffness matrix $Q_H^*$ along axes $L$ and $T$. The material reduced stiffness matrix $Q_H^*$ is defined with the five homogenized parameters of the honeycomb: $E_L^H$, $E_T^H$, $G_{LT}^H$, $\nu_{LT}^H$, $\nu_{TL}^H$.

$$ Q_I = \begin{pmatrix} \frac{E}{1-\nu^2} & \frac{\nu E}{1-\nu^2} & 0 \\ \frac{\nu E}{1-\nu^2} & \frac{E}{1-\nu^2} & 0 \\ 0 & 0 & G \end{pmatrix} \quad \text{and} \quad Q_H = Q_H^* = \begin{pmatrix} \frac{E_L^H}{1-\nu_{LT}^H \nu_{TL}^H} & \frac{\nu_{LT}^H E_T^H}{1-\nu_{LT}^H \nu_{TL}^H} & 0 \\ \frac{\nu_{LT}^H E_T^H}{1-\nu_{LT}^H \nu_{TL}^H} & \frac{E_L^H}{1-\nu_{LT}^H \nu_{TL}^H} & 0 \\ 0 & 0 & G_{LT}^H \end{pmatrix} \quad (6) $$

The bending matrix of the total sandwich is represented by $D^F$:

$$ D^F = \begin{pmatrix} D_{11} & D_{12} & D_{16} \\ D_{12} & D_{22} & D_{26} \\ D_{16} & D_{26} & D_{66} \end{pmatrix}. \quad (7) $$

All of the elements of the bending matrix $D^F$ are defined with the reference reduced stiffness matrices of each material $Q_I$ and $Q_H$. The square matrices $D^F$, $Q_I$ and $Q_H$ with nine terms can be replaced respectively by simple column $\{D_{\alpha\beta}\}$, $\{Q_{\alpha\beta}\}_I$ and $\{Q_{\alpha\beta}\}_H$ with the only six terms $\alpha\beta$: 11, 22, 12, 66, 16 and 26. With such a representation, the column $\{D_{\alpha\beta}\}$ can be expressed
with the columns {$Q_{\alpha\beta}$}$_{I}$ and {$Q_{\alpha\beta}$}$_{H}$ by

$$
\{D_{\alpha\beta}\}_{SHI} = \frac{2}{3}[d_1 \cdot \{Q_{\alpha\beta}\}_I + d_0 \cdot \{Q_{\alpha\beta}\}_H], \quad \text{where } d_0 = \frac{h_H^3}{8} \text{ and } d_1 = \frac{3h_H^2 e + 3h_H^3 e^2 + e^3}{4} \tag{8}
$$

SHI is the abbreviation for the symmetric honeycomb sandwich with isotropic faces.

$d_0$ is a term associated with the half-honeycomb central ply and $d_1$ is related to the isotropic face. The development of each term of {$D_{\alpha\beta}$} leads to the usual formulations:

$$
D_{11} = \frac{h_H^2}{2} e + h_H e^2 + \frac{2}{3} e^3 \left( \frac{E}{(1-\nu^2)} + \frac{h_H^3}{12} \frac{E_L^H}{(1-\nu_{LT}^H \nu_{TL}^H)} \right), \quad D_{22} = \frac{h_H^2}{2} e + h_H e^2 + \frac{2}{3} e^3 \left( \frac{E}{(1-\nu^2)} + \frac{h_H^3}{12} \frac{E_T^H}{(1-\nu_{LT}^H \nu_{TL}^H)} \right),
$$

$$
D_{12} = \frac{h_H^2}{2} e + h_H e^2 + \frac{2}{3} e^3 \left( \frac{\nu E}{(1-\nu^2)} + \frac{h_H^3}{12} \frac{\nu_{LT}^H E_T^H}{(1-\nu_{LT}^H \nu_{TL}^H)} \right), \quad D_{66} = \frac{h_H^2}{2} e + h_H e^2 + \frac{2}{3} e^3 \left( G + \frac{h_H^3}{12} G_{LT}^H \right), \quad D_{16} = 0, \quad D_{26} = 0. \quad (9)
$$

We can frequently neglect the honeycomb properties in comparison with the face properties. Therefore, if we discard the
second term of each element of the stiffness bending matrix, it leads to a simplification of all terms. Using this condition, the
sandwich has the same simple bending rigidity D in the total plane xy. Then, the sandwich has an equivalent isotropic
behavior stiffness:

$$
D = \frac{h_H^2}{2} e + h_H e^2 + \frac{2}{3} e^3 \left( \frac{E}{(1-\nu^2)} \right). \qquad (10)
$$

2.4. Transverse shear stiffness

For the reference transverse shear stiffness S, we can use Allen formula Eq. (2), but it does not work when the honeycomb thickness approaches zero. A few years ago, Gay [11] developed a model for the transverse shear moduli of a beam composed of three layers. The face thickness must be small compared to the core thickness. The materials are isotropic or transverse isotropic. We have adapted the model to a symmetric honeycomb plate with equal face sheets. We have established two equivalent transverse shear moduli $G_{xzeq}$ and $G_{xyeq}$, expressed as follows:

$$
G_{xzeq} = \frac{1}{k_{xz}} \cdot \frac{2eG + h_H G_{xz}^H}{(h_H + 2e)} \quad \text{and} \quad G_{yzeq} = \frac{1}{k_{yz}} \cdot \frac{2eG + h_H G_{yz}^H}{(h_H + 2e)}, \qquad (11)
$$

with the following Gay nomenclature

$$
\begin{align}
k_{xz} \text{ or } yz &= \frac{a}{8[E_2 H_2^3 + E_1 (H_1^3 - H_2^3)]} \left\{ \frac{E_2 H_2^3}{G_2} \left[ E_1 H_1^2 + \left(\frac{4}{5}E_2 - E_1\right) H_2^2 \right] + \frac{E_1^2}{G_1} \frac{4}{5} H_1^5 + \frac{H_2^5}{5} - H_1^2 H_2^3 \right\} \nonumber \\
&\quad + \frac{3bE_1 (H_1^2 - H_2^2)}{E_2 H_2^3 + E_1 (H_1^3 - H_2^3)} \nonumber \\
a &= \frac{\langle GS \rangle}{\langle EI \rangle} = 12 \frac{G_2 H_2 + G_1 (H_1 - H_2)}{E_2 H_2^3 + E_1 (H_1^3 - H_2^3)}, \quad b = \frac{a}{16} H_2 \frac{E_1}{G_1} \left\{ \frac{H_2^2}{3} + H_1^2 \left( \frac{G_1}{G_2} - 1 \right) - H_2^2 \frac{G_1}{G_2} \left( 1 - \frac{2E_2}{3E_1} \right) \right\} \nonumber \\
H_1 &= h_H + 2e, \quad H_2 = h_H, \quad G_1 = G_{\text{face}} = G = \frac{E}{2(1+\nu)}, \nonumber \\
G_2 &= G_{\text{Honeycomb}} = G_{xz}^H \text{ or } G_{yz}^H \quad \text{(following the analyzed reference plane xz or yz)} \tag{12}
\end{align}
$$

The material *L*, *T* and *T'* axes of the honeycomb are respectively parallel to the reference *x*, *y* and *z* axes of the plate. $G_{xz}^H$ and $G_{yz}^H$ of the honeycomb are thus simply the material moduli $G_{LT}^H$ and $G_{TT'}^H$.

The reference transverse shear stiffness *S* being homogenized on the multilayer, it is analytically expressed with *h*<sub>s</sub> and
the two equivalent transverse shear moduli:

$$
S = h_S (G_{xzeq} \cdot G_{yzeq})^{1/2}. \qquad (13)
$$

It is important to note that, with the conditions $h_H \gg 2e$, $G \gg G_{xz}^H$ and $G \gg G_{yz}^H$, the analytical $G_{xzeq}$ and $G_{xyeq}$ converge simply and respectively towards:

$$
G_{xzeq}^{\text{simpl}} \approx \frac{(h_H + e)^2}{h_H \cdot h_S} G_{xz}^H \quad \text{and} \quad G_{yzeq}^{\text{simpl}} \approx \frac{(h_H + e)^2}{h_H \cdot h_S} G_{yz}^H. \qquad (14)
$$
Using Eqs. (13) and (14), S can be written as (15). S is rigorously identical to the analytical Allen formula (2), with the honeycomb instead of the core:

$$S = h_H \left(1 + \frac{e}{h_H}\right)^2 G_g^H \quad \text{with} \quad G_g^H = (G_{xz}^H \cdot G_{yz}^H)^{1/2}. \qquad (15)$$

For the rest of the presentation, the terms $G_{xzeq}^{\text{simpl}}$ and $G_{yzeq}^{\text{simpl}}$ are renamed as $G_{\text{Allen}}$ and $G_{\text{Allen}}$.

## 2.5. Analytical wavenumber modulus

With the two Laplacian operators $\nabla^2$ and $\nabla^4$, the differential equation (3) is rewritten as

$$D\left(\frac{\partial^4 w}{\partial x^4} + 2\frac{\partial^4 w}{\partial x^2 \partial y^2} + \frac{\partial^4 w}{\partial y^4}\right) + m_s \frac{\partial^2 w}{\partial t^2} - m_s \frac{D}{S} \frac{\partial^2}{\partial t^2} \left(\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2}\right) = 0. \qquad (16)$$

Using the relation between the acceleration and the displacement (for harmonic motions), we obtain

$$D\left(\frac{\partial^4 w}{\partial x^4} + 2\frac{\partial^4 w}{\partial x^2 \partial y^2} + \frac{\partial^4 w}{\partial y^4}\right) - m_s \omega^2 w + m_s \frac{D}{S} \omega^2 \left(\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2}\right) = 0. \qquad (17)$$

For an infinite plate, the expression of a propagating wave is $\tilde{W}(x, y, t) = w_0 \cdot e^{j(\omega t - k_x x - k_y y)}$, where $k_x$ and $k_y$ are the wavenumber components in the x and y directions. The wavenumber modulus is related to $k_x$ and $k_y$ by the relation $k^2 = k_x^2 + k_y^2$. We derive the relations of $\nabla^2$ and $\nabla^4$ with the wavenumber in the directions x and y:

$$\nabla_{\text{plate}}^2 w = -k^2 w, \quad \nabla_{\text{plate}}^4 w = +k^4 w. \qquad (18)$$

With the differential equation, the wavenumber equation is finally

$$\left[ D(k_x^2 + k_y^2)^2 - m_s \omega^2 - m_s \frac{D}{S} (k_x^2 + k_y^2) \omega^2 \right] w = 0 \quad \text{or} \quad Dk^4 - m_s \omega^2 - m_s \frac{D}{S} k^2 \omega^2 = 0. \qquad (19)$$

The modulus of the wavenumber $k(\omega)$ is extracted from Eq. (19). It is a function of ($m_s$, $D$, $S$, $\omega$).

$$|k(\omega)| = \left( \left( \frac{1}{2S} \right) \left( m_s \omega^2 + \omega \left( m_s^2 \omega^2 + \left( 4m_s^2 S^2 / D \right) \right) \right)^{1/2} \right)^{1/2} \qquad (20)$$

## 2.6. Analytical modal frequencies under simply supported boundary conditions

The circular frequency $\omega$ can also be expressed from Eq. (19) as

$$\omega = \left(\frac{D}{m_s}\right)^{1/2} (k_x^2 + k_y^2) \left(1 + \frac{D}{S}(k_x^2 + k_y^2)\right)^{-1/2}. \qquad (21)$$

For a rectangular plate, under simply supported boundary conditions along the four sides and having the dimensions ($a,b$), $k_x$ and $k_y$ can be expressed with the modal numbers $i$ and $j$ associated to the numbers of a half-stationary sine wave:

$$k_x = \frac{i\pi}{a}, \quad k_y = \frac{j\pi}{b}. \qquad (22)$$

The solution for the bending displacement may be written as

$$W_{ij}(x, y) = w_{0 \ i j} \sin(k_x x) \sin(k_y y) = w_{0 \ i j} \sin\left(\frac{i\pi}{a}x\right) \sin\left(\frac{j\pi}{b}y\right).$$

Introducing Eq. (22) into Eq. (21), and knowing $\omega=2\pi f$, the modal frequencies of a simply supported rectangular sandwich plate with transverse shear are

$$f_{\text{SSP},i,j}^{\text{SHI}} = \frac{1}{2\pi} \left( \frac{D}{m_s} \right)^{1/2} \left( \left( \frac{i\pi}{a} \right)^2 + \left( \frac{j\pi}{b} \right)^2 \right) \left( 1 + \frac{D}{S} \left( \left( \frac{i\pi}{a} \right)^2 + \left( \frac{j\pi}{b} \right)^2 \right) \right)^{-1/2}. \qquad (23)$$

SSP is the abbreviation for the plate under simply supported boundary conditions along its four sides. In this case the boundary conditions are

- along the two edges a: the displacements $u_x=v_y=w_z=0$ and the rotations $R_y=R_z=0$,

- along the two edges b: the displacements $u_x=v_y=w_z=0$ and the rotations $R_x=R_z=0$.

If $h_N$ equals zero the sandwich becomes a simple plate with only two isotropic face sheets. In this case, S is classically much higher than D. Therefore, Eq. (23) becomes again the isotropic plate formula.
## 2.7. Modal density for all boundary conditions

From Eq. (23), we can easily determine the modal density using the average frequency interval between two successive modes $\bar{\delta}(f)$, taken over a frequency bandwidth:

$$n(f) = \frac{1}{\bar{\delta}(f)}. \qquad (24)$$

However, this would be valid only for a simply supported plate. In statistical energy analysis (SEA), the modal density is defined from the wavenumber, through the phase and group celerities of bending waves in the plate:

$$ \left. \begin{array}{l} \text{Phase celerity: } c_p = \frac{\omega}{k} \\ \text{Group celerity: } c_g = \frac{\partial\omega}{\partial k} \\ \text{Area of the plate: } A \end{array} \right\} \bar{\delta}(\omega) = \frac{c_p c_g}{\omega A} \Rightarrow n(\omega) = \frac{1}{\bar{\delta}(\omega)}. \qquad (25) $$

Clarkson and Ranky [5] used this principle and yielded the modal density for a honeycomb sandwich panel, independently from any boundary conditions:

$$ n(\omega)^{SHI}_{\text{boundary conditions}} = \frac{A m_s \omega}{2S} \left( 1 + m_s \omega^2 + \frac{2S^2}{D} \right) m_s^2 \omega^4 + \frac{4m_s \omega^2 S^2}{D} \left( -1/2 \right). \qquad (26) $$

Like the wavenumber modulus, the modal density is also a function of the same four parameters: the mass per unit area $m_s$, the bending rigidity $D$, the out-of-plane shear stiffness $S$ and the circular frequency $\omega$.

## 2.8. Mechanical impedance

The modal density is currently obtained by experimental techniques, with the spatial average of the real part of the frequency mobility. The mobility of a structure, also called the admittance, is the ratio of the dynamic velocity to the force, at the excitation point. The modal density is written as

$$ n(\omega) = 4m_s A \left\langle \text{Re} \left( \frac{V}{F}(\omega) \right) \right\rangle^{\text{at excitation points}}_{\text{finite plate}}. \qquad (27) $$

Term A represents the area of the finite plate. The external brackets indicate the spatial average. The spatial average of the real part of the frequency response function $(V/F)$ is also considered equivalent to the inverse mechanical impedance of the infinite plate. Introducing the modal density from Eq. (26) into Eq. (27), therefore an analytical expression of the inverse infinite mechanical impedance is obtained as

$$ \begin{aligned} Z_{\infty}^{-1}(\omega) &= \left\langle \text{Re} \left( \frac{V}{F}(\omega) \right) \right\rangle^{\text{at excitation point s}}_{\text{finite plate}} \\ &= \frac{\omega}{8S} \left( 1 + m_s \omega^2 + \frac{2S^2}{D} \right) m_s^2 \omega^4 + \frac{4m_s \omega^2 S^2}{D} \left( -1/2 \right). \end{aligned} \qquad (28) $$

## 2.9. Transition frequency

When $m_s$, $D$, and $S$ are fixed, the wavenumber module, modal density and mechanical impedance functions are only a function of the circular frequency $\omega$. A transition frequency appears and separates all curves into two frequency domains:

$$ f_{\text{trans}} = \frac{S}{4\pi} \left( \frac{1}{m_s D} \right)^{1/2}. \qquad (29) $$

In the first domain, below the transition frequency (in the low frequency domain), the orthotropic sandwich panel has a classical isotropic plate behavior. In the second domain, above the transition frequency (in the high frequency domain), the out-of-plane shear rigidity is very significant and changes the behavior. Table 1 gives simplified formulations for the two frequency domains.

## 2.10. Critical frequency

The critical frequency is important for the acoustic behavior of a plate. Under a given mechanical excitation, the sound radiated from a finite plate is maximum at this critical frequency and then decreases to a limit value for higher frequencies. Under a diffused acoustic excitation, the panel transmission loss is particularly weak at this same frequency, which is the lowest of the coincidence frequencies. The coincidence frequencies change with the angle of the incident acoustic wave. The
<table><caption>**Table 1**

Simplified formulations for the two frequency domains.</caption><thead><tr><th></th><th colspan="2" style="text-align:center;">f<sub>trans</sub> = <sup>S</sup>&frasl;<sub>4&zeta;</sub>(<sup>m</sup>&frasl;<sub>D</sub>)<sup>1/2</sup></th></tr><tr><th></th><th style="text-align:center;">Below transition frequency</th><th style="text-align:center;">Above transition frequency</th></tr></thead><tbody><tr><td>Wavenumber modulus k</td><td>k(f)<sub>f→0</sub> = (<sup>m</sup>&frasl;<sub>D</sub>)<sup>1/4</sup>(2&pi;f)<sup>1/2</sup></td><td>k(f)<sub>f→&infin;</sub> = (<sup>m</sup>&frasl;<sub>S</sub>)<sup>1/2</sup>2&pi;f</td></tr><tr><td>Phase celerity c<sub>p</sub> = <sup>&omega;</sup>&frasl;<sub>k</sub></td><td>c<sub>p</sub>(f)<sub>f→0</sub> = (<sup>D</sup>&frasl;<sub>m</sub>)<sup>1/4</sup>(2&pi;f)<sup>1/2</sup></td><td>c<sub>p</sub>(f)<sub>f→&infin;</sub> = (<sup>S</sup>&frasl;<sub>m</sub>)<sup>1/2</sup></td></tr><tr><td>Group celerity c<sub>g</sub> = <sup>d&alpha;</sup>&frasl;<sub>dk</sub></td><td>as c<sub>p</sub> = cte&middot;&omega;<sup>1/2</sup> with cte = (D/m<sub>s</sub>)<sup>1/4</sup><br>&rarr; c<sub>g</sub> = 2c<sub>p</sub></td><td>as c<sub>p</sub> = cte with cte = (S/m<sub>s</sub>)<sup>1/2</sup><br>&rarr; c<sub>g</sub> = c<sub>p</sub></td></tr><tr><td>Modal density n<sub>plate</sub> = <sup>d&alpha;</sup>&frasl;<sub>c<sub>p</sub>c<sub>g</sub></sub></td><td>n(f)<sub>f→0</sub> = <sup>4</sup>&frasl;<sub>2</sub>(<sup>m</sup>&frasl;<sub>D</sub>)<sup>1/2</sup></td><td>n(f)<sub>f→&infin;</sub> = A(<sup>m</sup>&frasl;<sub>S</sub>)&middot;2&pi;f</td></tr><tr><td>Spatial real mobility V/F (admittance) &lang;Re(V/F)&rang; = <sup>n<sub>plate</sub></sup>&frasl;<sub>4m&middot;A</sub></td><td>&lang;Re(V/f)<sub>f→0</sub>&rang; = <sup>1</sup>&frasl;<sub>8(Dm<sub>s</sub>)<sup>1/2</sup></sub></td><td>&lang;Re(V(f)<sub>f→&infin;</sub>&rang; = (<sup>&alpha;</sup>&frasl;<sub>S</sub>)f</td></tr><tr><td>Mechanical impedance Z<sub>&infin;</sub> = 1/&lang;Re(V/F)&rang;</td><td>Z<sub>&infin;</sub>(f)<sub>f→0</sub> = 8(Dm<sub>s</sub>)<sup>1/2</sup></td><td>Z<sub>&infin;</sub>(f)<sub>f→&infin;</sub> = <sup>2S</sup>&frasl;<sub>Sf</sub></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Symmetric honeycomb sandwich with composite faces.</figcaption></figure>

lowest of the coincidence frequencies and critical frequency are linked. This frequency is located at the intersection of the acoustic wavenumber function with the bending wavenumber function of the material that the panel is composed of. The bending wavenumber function for the homogenized sandwich is represented by Eq. (20). The acoustic wavenumber is simply defined by

$$k_{ac}(\omega) = \frac{\omega}{c_{air}}. \qquad (30)$$

The intersection of these two functions gives the critical frequency

$$f_{cr} = \left(\frac{m_s}{D}\right)^{1/2} \frac{c_{air}^2}{2\pi} \left(1 - \frac{m_s c_{air}^2}{S}\right)^{-1/2}. \qquad (31)$$

In the critical frequency formulation, $m_s$, $D$ and $S$ can be seen again.

### 3. Symmetric honeycomb sandwich with composite faces

We are now going to analyze a symmetric honeycomb sandwich with composite faces, shown in Fig. 2, and to compare it with the symmetric honeycomb sandwich with isotropic faces in Fig. 1.

A few years after Clarkson and Ranky, authors Renji et al. [6] worked on the modal densities of honeycomb sandwich panels with composite faces. The honeycomb core and the composite faces are both orthotropic. To include the effect of transverse shear, Mindlin's theory is used. In this new case, the classical laminate theory (CLT) is necessary to determine the membrane, bending and transverse rigidity matrices of the total laminate. In order to eliminate the eventual coupling between membrane and flexion, the sandwich is also assumed to be symmetric here.

The stiffness and inertia parameters may be established by considering the complete plate as a multilayer. From the bottom to the top (see Fig. 3), the plies are designated as layers 1 to n. The total thickness of the sandwich is $h_s$. The central layer represents the honeycomb core with thickness $h_H$. Each layer identified by index k is bounded in the thickness direction z by the coordinates $z_{k-1}$ and $z_k$. Each ply is characterized by the thickness $e_k$, the density $\rho_k$, the reduced stiffness matrix $Q_{\alpha\beta}^{*k}$ and the two transverse shear moduli expressed in the material axes, the orientation $\theta_k$ of the material axes relative to reference axes.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Laminate.</figcaption></figure>

### 3.1. Differential equation of motion

Some authors, such as Lesueur [15] and, more recently, Yu and Cleghorn [9], suggest that the first $I_1$ and third $I_3$ order mass moments of inertia be introduced in the differential equation of motion. We are able to generalize the differential equation (16) and that proposed by Renji and Nair for an orthotropic plate with transverse shear [6], further including $I_1$ and $I_3$, as the following equation:

$$ \left( D_{11} \frac{\partial^4 w}{\partial x^4} + 2(D_{12} + 2D_{66}) \frac{\partial^4 w}{\partial x^2 \partial y^2} + D_{22} \frac{\partial^4 w}{\partial y^4} \right) + I_1 \frac{\partial^2 w}{\partial t^2} - I_3 \frac{\partial^2}{\partial t^2} \left( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \right) - \frac{I_1}{S} \frac{\partial^2}{\partial t^2} \left( D_{11} \frac{\partial^2 w}{\partial x^2} + D_{22} \frac{\partial^2 w}{\partial y^2} \right) = 0. \quad (32) $$

### 3.2. Mass inertia moments $I_1$ and $I_3$

$I_1$ and $I_3$ represent the first- and third-order mass moments of inertia for the laminate:

$$ I_1 = \sum_{k=1}^{n} \int_{z_{k-1}}^{z_k} \rho_k \, dz = \sum_{k=1}^{n} \rho_k (z_k - z_{k-1}) = \sum_{k=1}^{n} \rho_k e_k, \quad I_3 = \sum_{k=1}^{n} \int_{z_{k-1}}^{z_k} \rho_k z^2 \, dz = \frac{1}{3} \sum_{k=1}^{n} \rho_k (z_k^3 - z_{k-1}^3). \quad (33) $$

$I_1$ is identical to the previous term $m_5$.

### 3.3. Bending stiffness

$D_{11}$, $D_{22}$, $D_{12}$, and $D_{66}$ are the four elements of the bending rigidity matrix $D^F$ of the total laminate. $D^F$ is a 3 × 3 matrix defined by Eq. (7). It is understood that $D_{16}$ and $D_{26}$ are assumed to have been neglected or to be equal to zero in Eq. (32). All of the individual elements $D_{\alpha\beta}$ are written as follows:

$$ D_{\alpha\beta} = \int_{-h_s/2}^{+h_s/2} z^2 Q_{\alpha\beta} \, dz = \sum_{k=1}^{n} \int_{z_{k-1}}^{z_k} z^2 Q_{\alpha\beta}^k \, dz = \frac{1}{3} \sum_{k=1}^{n} (z_k^3 - z_{k-1}^3) Q_{\alpha\beta}^k. \quad (34) $$

$Q_{\alpha\beta}^k$ are the elements of the reference reduced stiffness matrix $Q_k$ of ply $k$, along the reference axes x and y. Reference reduced stiffness matrix $Q_k$ of ply $k$ along reference axes

$$ Q_k = \begin{pmatrix} Q_{11}^k & Q_{12}^k & Q_{16}^k \\ Q_{12}^k & Q_{22}^k & Q_{26}^k \\ Q_{16}^k & Q_{26}^k & Q_{66}^k \end{pmatrix} \quad (35) $$

They are calculated by taking into account the elements of the material reduced stiffness matrix of ply $k$ along material axes $Q_{\alpha\beta}^k$ and the rotation angle $\theta_k$ of the ply $k$.

Material reduced stiffness matrix $Q_k^*$ of ply $k$ along material axes

$$ Q_k^* = \begin{pmatrix} Q_{11}^{*k} & Q_{12}^{*k} & 0 \\ Q_{12}^{*k} & Q_{22}^{*k} & 0 \\ 0 & 0 & Q_{66}^{*k} \end{pmatrix} \quad (36) $$
By applying the transformation **T** matrix simply to the elements of the matrices $Q_k^*$ presented in column form $\{Q_{\alpha\beta}^*\}_k$, the following process can be used:

$$ \{Q_{\alpha\beta}\}_k = \mathbf{T} \cdot \{Q_{\alpha\beta}^*\}_k, \qquad (37) $$

with

$$ \begin{pmatrix} Q_{11}^k \\ Q_{22}^k \\ Q_{12}^k \\ Q_{66}^k \\ Q_{16}^k \\ Q_{26}^k \end{pmatrix} = \begin{pmatrix} c^4 & s^4 & 2c^2s^2 & 4c^2s^2 \\ s^4 & c^4 & 2c^2s^2 & 4c^2s^2 \\ c^2s^2 & c^2s^2 & (c^4+s^4) & -4c^2s^2 \\ c^2s^2 & c^2s^2 & -2c^2s^2 & (c^2-s^2)^2 \\ c^3s & -cs^3 & cs(s^2-c^2) & 2cs(s^2-c^2) \\ cs^3 & -c^3s & cs(c^2-s^2) & 2cs(c^2-s^2) \end{pmatrix} \cdot \begin{pmatrix} Q_{11}^{*k} \\ Q_{22}^{*k} \\ Q_{12}^{*k} \\ Q_{66}^{*k} \end{pmatrix} \qquad (38) $$

and where $c=\cos \theta_k$ and $s=\sin \theta_k$.

### 3.4. Transverse shear stiffness

$S$ is the global transverse shear rigidity of the laminate and now depends on the two homogenized transverse shear elements called $C_{11}^{\text{hom}}$ and $C_{22}^{\text{hom}}$ of the $\mathbf{C}^{\text{hom}}$ matrix. They replace $G_{xzeq}$ and $G_{xyeq}$ into Eq. (13) to give

$$ S = h_s (C_{11}^{\text{hom}} \cdot C_{22}^{\text{hom}})^{1/2}. \qquad (39) $$

Berthelot [12] has developed formulations for $C_{11}^{\text{hom}}$, $C_{22}^{\text{hom}}$ and $C_{12}^{\text{hom}}$ of a laminate. The symmetric homogenized transverse shear $\mathbf{C}^{\text{hom}}$ is

$$ \mathbf{C}^{\text{hom}} = \begin{pmatrix} C_{11}^{\text{hom}} & C_{12}^{\text{hom}} \\ C_{12}^{\text{hom}} & C_{22}^{\text{hom}} \end{pmatrix}, \qquad (40) $$

and has in fact three elements:

$$ C_{11}^{\text{hom}} = k_{11} \frac{\sum_{k=1}^{n} (G_{23}^{k} \sin^2 \theta_k + G_{13}^{k} \cos^2 \theta_k) e_k}{h_s}, \qquad (41) $$

$$ C_{22}^{\text{hom}} = k_{22} \frac{\sum_{k=1}^{n} (G_{23}^{k} \cos^2 \theta_k + G_{13}^{k} \sin^2 \theta_k) e_k}{h_s}, \qquad (42) $$

$$ C_{12}^{\text{hom}} = k_{12} \frac{\sum_{k=1}^{n} (G_{13}^{k} - G_{23}^{k}) \cos \theta_k \sin \theta_k e_k}{h_s}, \qquad (43) $$

$k_{11}$, $k_{22}$ and $k_{12}$ are correction factors of the transverse shear of the total laminate. If a ply $k$ is composed of isotropic or orthotropic material, with material axes parallel to the reference axes ($\theta_k=0^\circ$), Table 2 shows how we obtain in planes $xz$

<table><caption>**Table 2**

Transverse shear moduli as a function of the nature of the ply.</caption><thead><tr><th colspan="2"></th><th colspan="2">Properties along the material axes</th></tr><tr><th colspan="2"></th><th colspan="2">Type of ply</th></tr><tr><th colspan="2">Isotropic</th><th colspan="2">Orthotropic</th></tr><tr><th colspan="2"></th><th>Unidirectional composite UD<br>Material axes (1,2,3)</th><th>Honeycomb<br>Material axes (L,T,T)</th></tr></thead><tbody><tr><td colspan="2">$G_{13}^k = G_{23}^k = G = E/2(1+\nu)$</td><td>$G_{13}^k = G_{13}^{UD} = G_{23}^{UD} = G_{23}^{UD}$<br>If the material axes are parallel to the<br>reference axes<br>$\theta_k = 0^\circ$</td><td>$G_{13}^k = G_{13}^{H} = G_{23}^{H} = G_{23}^{H}$<br>If the material axes are parallel to the<br>reference axes<br>$\theta_k = 0^\circ$</td></tr><tr><td><b>Plane</b></td><td><b>Formulations</b></td><td colspan="2"><b>Properties along the reference axes (x,y,z)</b></td></tr><tr><td><b>xz</b></td><td>$G_{23}^k \sin^2 \theta_k + G_{13}^k \cos^2 \theta_k$</td><td>G &forall; $\theta_k$</td><td>$G_{13}^k$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $G_{13}^{UD}$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $G_{13}^{H}$</td></tr><tr><td><b>yz</b></td><td>$G_{23}^k \cos^2 \theta_k + G_{13}^k \sin^2 \theta_k$</td><td>G &forall; $\theta_k$</td><td>$G_{23}^k$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $G_{23}^{UD}$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $G_{23}^{H}$</td></tr><tr><td colspan="2">$(G_{13}^k - G_{23}^k) \cos \theta_k \sin \theta_k$</td><td>0 &forall; $\theta_k$</td><td>$0$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $0$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $0$</td></tr></tbody></table>
and yz, a classical shear for isotropic material and the material transverse shear properties for an unidirectional composite (UD) and a honeycomb.

For a multilayer with only isotropic plies of the same material, we obtain $C_{11}^{\text{hom}} = k_{11}G$, $C_{22}^{\text{hom}} = k_{22}G$ and $C_{12}^{\text{hom}} = 0$. The transverse shear rigidity is equal to $S = h_S(k_{11} \cdot k_{22})^{1/2}G$. In this case, $k_{11}$ and $k_{22}$ have the same transverse shear coefficient $k_c$. The shear coefficient is used to account for the variation of the shear stress across the cross section. It represents a kind of ratio of the mean stress to the maximum stress. Its value depends on the beam or plate (with rectangular cross section) theory used by Timoshenko (2/3), Mindlin ($\pi^2/12$), Reissner (5/6), Cowper, Stephen (see [16,17]), Hutchinson (see [18]), etc. Finally, for an isotropic plate, we find the general relation $S = h_S k_c G$. For a symmetric multilayer with honeycomb core and composite faces, the terms of the transverse shear matrix are more complicated but sometimes they can be simplified as a sandwich with isotropic faces. This point will be discussed later, with an application, in Section 4.3.

Like we have done for a symmetric honeycomb sandwich with isotropic faces, we are able to analyze the wavenumber, the modal frequencies and the modal density of a symmetric honeycomb sandwich with composite faces.

### 3.5. Wavenumber modulus

From the differential equation (32) and from Eq. (18), the wavenumber is a solution of

$$ (D_{11}k_x^4 + 2(D_{12} + 2D_{66})k_x^2k_y^2 + D_{22}k_y^4) - I_1\omega^2 - I_3(k_x^2 + k_y^2)\omega^2 - I_1/S(D_{11}k_x^2 + D_{22}k_y^2)\omega^2 = 0. \quad (44) $$

### 3.6. Analytical modal frequencies under simply supported conditions

An expression of the circular frequency $\omega$ is extracted From Eq. (44)

$$ \omega = (D_{11}k_x^4 + 2(D_{12} + 2D_{66})k_x^2k_y^2 + D_{22}k_y^4)^{1/2} (I_1 + I_3(k_x^2 + k_y^2) + I_1/S(D_{11}k_x^2 + D_{22}k_y^2))^{-1/2}. \quad (45) $$

For a rectangular plate, simply supported and having the dimensions (a,b), $k_x$ and $k_y$ can be expressed with the modal numbers i and j associated to the half-stationary sine wavenumber, Eq. (22). The modal frequencies of a rectangular symmetric honeycomb sandwich plate with composite faces and simply supported with transverse shear are

$$ f_{\text{SSP},i,j}^{\text{SHC}} = \frac{1}{2\pi} \left(\frac{1}{I_1}\right)^{1/2} \left( \frac{D_{11}(i\pi/a)^4 + 2H(i\pi/a)^2(j\pi/b)^2 + D_{22}(j\pi/b)^4}{1 + \frac{I_3}{I_1}\left((i\pi/a)^2 + (j\pi/b)^2\right) + \frac{1}{S}\left(D_{11}(i\pi/a)^2 + D_{22}(j\pi/b)^2\right)} \right)^{1/2}, \quad (46) $$

with

$$ H = (D_{12} + 2D_{66}). \quad (47) $$

SHC is the abbreviation for the symmetric honeycomb sandwich with composite faces. If the composite faces are a multilayer made up of the same unidirectional ply material and thickness, we obtain the two mass moments of inertia of the laminate:

$$ I_1 = \rho_H h_H + 2\rho_{UD} e = m_s, \quad I_3 = \frac{\rho_H h_H^3}{12} + \frac{\rho_{UD}(h_S^3 - h_H^3)}{12}. \quad (48) $$

Eq. (46) is mostly used under two kinds of condition examined just after.

### 3.7. First condition: S very high

The first condition is when $I_3/I_1 \gg D_{11}/S$ and $D_{22}/S$. This appears when the two transverse shear moduli of the total laminate are so high that the shear rigidity $S$ is very high in front of the bending stiffness terms. The effect of the transverse shear can be neglected. This condition is usually used for the classical plate theory (CPT), which applies to a thin orthotropic plate without transverse shear. The modal frequencies formula simply becomes

$$ f_{\text{SSP},i,j}^{\text{CPT}} = \frac{1}{2\pi} \left(\frac{1}{I_1}\right)^{1/2} \left( D_{11} \left(\frac{i\pi}{a}\right)^4 + 2H \left(\frac{i\pi}{a}\right)^2 \left(\frac{j\pi}{b}\right)^2 + D_{22} \left(\frac{j\pi}{b}\right)^4 \right)^{1/2} \left( 1 + \frac{I_3}{I_1} \left( \left(\frac{i\pi}{a}\right)^2 + \left(\frac{j\pi}{b}\right)^2 \right) \right)^{-1/2}. \quad (49) $$
### 3.8. Second condition: S very low

#### 3.8.1. Simplified analytical modal frequencies under simply supported conditions

Opposite from the first condition, we find the case where $I_3/I_1 \ll D_{11}/S$ and $D_{22}/S$. This condition occurs when a honeycomb is present in the core. The two transverse shear moduli of the total laminate are so low that the transverse shear stiffness S is very low in front of the bending stiffness terms. The effect of transverse shear is more significant than the two mass moments of inertia of the laminate. The modal frequencies formulation is now

$$f_{\text{SSP},i,j}^{\text{SHC with N low}} = \frac{1}{2\pi} \left(\frac{1}{I_1}\right)^{1/2} \left( D_{11} \left(\frac{i\pi}{a}\right)^4 + 2H \left(\frac{i\pi}{a}\right)^2 \left(\frac{j\pi}{b}\right)^2 + D_{22} \left(\frac{j\pi}{b}\right)^4 \right)^{1/2} \left(1 + \frac{1}{S} \left(D_{11} \left(\frac{i\pi}{a}\right)^2 + D_{22} \left(\frac{j\pi}{b}\right)^2\right)\right)^{-1/2} . \quad (50)$$

#### 3.8.2. Simplified analytical wavenumber modulus

Always under this last condition, the wavenumber equation (44) is reduced to

$$(D_{11}k_x^4 + 2(D_{12} + 2D_{66})k_x^2k_y^2 + D_{22}k_y^4) - I_1\omega^2 - I_1/S(D_{11}k_x^2 + D_{22}k_y^2)\omega^2 = 0. \quad (51)$$

#### 3.8.3. Simplified analytical modal density

The wavenumber equation above has been studied by Renji et al. [6], where $I_1$ is replaced by $m_s$. From this and from the SEA procedure, the authors have derived a general solution for the modal density $n(\omega)$. This time with the introduction of $D_{11}$, $D_{22}$, H (in-plane rigidities) and the out-of-plane shear rigidity S, the modal density is written as

$$n(\omega)^{\text{SHC with N low}} = \frac{Am_s\omega}{\pi S} \int_0^{\pi/2} \left\{ \frac{f_{2(\theta)}}{f_{1(\theta)}} + \frac{1}{f_{1(\theta)}} m_s\omega^2 f_{2(\theta)}^2 + \frac{2S^2}{(D_{11}D_{22})^{1/2}f_{1(\theta)}} \right\} m_s^2\omega^4 f_{2(\theta)}^2 + \frac{4m_s\omega^2S^2}{(D_{11}D_{22})^{1/2}f_{1(\theta)}} \right\} d\theta \quad (52)$$

where

$$\left\{ \begin{aligned} \gamma_1^2 &= \frac{1}{2} \left[ 1 - \frac{(D_{12} + 2D_{66})}{(D_{11}D_{22})^{1/2}} \right] \\ f_{1(\theta)} &= 1 - \gamma_1^2 \sin(2\theta) \\ f_{2(\theta)} &= \left( \frac{D_{11}}{D_{22}} \right)^{1/4} \cos^2(\theta) + \left( \frac{D_{22}}{D_{11}} \right)^{1/4} \sin^2(\theta) \end{aligned} \right. . \quad (53)$$

#### 3.8.4. Second condition, very low S: with $D_{11}=D_{22}=D$

##### 3.8.4.1. Simplified analytical wavenumber modulus. Renji et al. [19] established a solution for this wavenumber module if $D_{11}=D_{22}=D$:

$$|k(\omega)| = \left( \left( \frac{1}{2S} \right) \left( \frac{m_s \omega^2 + \omega \left( m_s^2 \omega^2 + \left( \frac{4m_s S^2}{D} \right) \left\{ \frac{3+(H/D)}{4} \right\} \right)^{1/2}}{\left\{ \frac{3+(H/D)}{4} \right\}} \right) \right)^{1/2} . \quad (54)$$

The ratio H/D represents the orthotropic coefficient. Now, the wavenumber modulus and the modal density depend on seven parameters: the mass per unit area $m_s$, the four terms of the bending rigidity matrix ($D_{11}$, $D_{22}$, $D_{12}$, $D_{66}$) through D and H, the transverse shear rigidity S and the circular frequency $\omega$.

#### 3.8.5. Second condition, very low S: with $D_{11}=D_{22}=D=H$

##### 3.8.5.1. Simplified analytical wavenumber modulus and modal density evolution. Using condition $D_{11}=D_{22}=D=H$, the ratio H/D becomes equal to unity. It is easy to see that the wavenumber modulus equation (54) proposed by Renji and Nair is the same as Eq. (20) proposed by Clarkson and Ranky. Still with the same conditions, the above Eq. (53) is simplified: $\gamma_1^2 = 0$, $f_{1(\theta)} = 1$, $\int_0^{\pi/2} d\theta = \pi/2$. The modal density of Eq. (52) becomes of course similar to Eq. (26), used for a symmetric honeycomb sandwich with isotropic faces.

Some particular arrangements in the plies can produce this situation. The initial complex symmetric honeycomb sandwich with composite faces returns to the same global properties ($m_s$, D, S) of the symmetric honeycomb sandwich with isotropic faces. Even though such sandwich seems very orthotropic at the beginning, with the presence of honeycomb and UD plies, in the end it can have isotropic plate rigidity. This can be analyzed later through the known formulations of the rigidity matrix. The difference of these two types of honeycomb sandwich (with isotropic or composites faces) in comparison with a real isotropic plate is the presence of the transverse shear rigidity. Table 1 shows above the transition
<figure><img src="image_6.png" /><figcaption>Fig. 4. Application of the symmetric honeycomb sandwich with composite faces.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Unit mesh composition of a honeycomb sandwich.</figcaption></figure>

frequency, the simultaneous influence of shear on wavenumber modulus, modal density and mechanical impedance. We will now analyze all of those points with our own application.

## 4. Application

In aerospace structures, most sandwich panels have light honeycomb and face sheets made of a fibrous composite panel. The application that has been developed is a symmetric honeycomb sandwich with faces as shown in Fig. 4. The central layer is a honeycomb made up of Nomex® material. The core is 12 mm thick. The skins are composed of eight plies of unidirectional composite. Each skin is a symmetric multilayer with angle arrangement [0°, 45°, 90°, -45°]s. The stack of one face is quasi-isotropic composite. The global skin thickness is 1 mm. Thus, the total sandwich is 14 mm thick. Two panel dimensions are used: panel A, 1.96(L) × 1.36(T) m² which is a rectangular plate and panel B, 1(L) × 1(T) m² which is a square plate.

### 4.1. Description of the properties

#### 4.1.1. Honeycomb

The Nomex® material of the honeycomb is paper composed of aramid fiber with phenolic resin. The basic properties of the web Nomex® material are as follows: density, 1240 kg/m³; Young modulus, 5.5 GPa; Poisson ratio, 0.33; Shear modulus, 2.07 GPa; thickness, 3 × 10⁻³ in. (76.2 µm).
The diameter of the hexagonal cell of the honeycomb $\phi_{cell}$ is 3/16 in. (4.76 mm). The material L, T and T' directions of the honeycomb ply are respectively parallel to the reference x, y and z directions of the panel. This means that there are double walls in the x direction and a single wall in the y direction (see Fig. 5). For the multilayer, this layer has $\theta_k=0^\circ$.

With these properties and using Gibson and Ashby formulations in [20] for a honeycomb with single and double walls, the material analytical properties derived for the honeycomb are

density $\rho_H = 53 \text{ kg/m}^3$ and engineering constants:

$$E_L = E_T = 0.27 \text{ MPa}, \quad \nu_{LT} = \nu_{TL} = 1, \quad G_{LT} = 0.162 \text{ MPa}, \quad E_T = 0.235 \text{ GPa},$$

$$\nu_{LT'} = \nu_{TT'} = 0, \quad \nu_{TL} = \nu_{TT} = 0.33, \quad G_{LT'} = 50.2 \text{ MPa}, \quad G_{TT'} = 33.1 \text{ MPa}.$$

The honeycomb layer is a structure, but is considered as a (fictious) homogeneous material. That is why some Poisson ratios may seem strange: $\nu_{LT}$ and $\nu_{TL}$ are theoretically equal to unity. However, for the computation of the material stiffness matrix, they are assumed to be close to 0.9, otherwise Eq. (9) cannot be computed. Since $E_T$ is large compared to $E_L$ and $E_T$, $\nu_{LT}$ and $\nu_{TT'}$ tend to zero. $\nu_{TL}$ and $\nu_{TT'}$ are equal to the web modulus.

### 4.1.2. Unidirectional composite (UD)

The unidirectional composite is an epoxy resin with high tenacity carbon yarn. The UD layer has the following basic properties along the material axes:

$$E_1 = 133.6 \text{ GPa}, \quad E_2 = 7.7 \text{ GPa}, \quad \nu_{12} = 0.29, \quad G_{12} = 3.1 \text{ GPa}, \\ \text{thickness } e_{UD} = 125 \text{ mm}, \quad \text{density } \rho_{UD} = 1550 \text{ kg/m}^3, \quad \text{fiber volumic ratio} = 55\%.$$

We consider the ply of UD as transverse isotropic, so the complementary properties are

$$E_3 = E_2 = 7.7 \text{ GPa}, \quad \nu_{13} = \nu_{12} = 0.29, \quad G_{13} = G_{12} = 3.1 \text{ GPa}, \\ \text{thus, when } \nu_{23} \text{ becomes } 0.5, \quad G_{23} = E_2/2(1+\nu_{23}) = 2.6 \text{ GPa}.$$

The other Young moduli are deduced from the reciprocity:

$$\nu_{21} = (E_2/E_1)\nu_{12}, \quad \nu_{31} = (E_3/E_1)\nu_{13} \quad \text{and} \quad \nu_{32} = (E_3/E_2)\nu_{23}.$$

The ply with angle zero is chosen in such a way that its material 1 and 2 directions are respectively in the reference x and y directions of the panel.

### 4.1.3. Material reduced stiffness matrices

The material reduced stiffness matrices $Q_H^*$ of the honeycomb ply is the same as Eq. (6) and $Q_{UD}^*$ for UD orthotropic plies is the following:

$$Q_{UD}^* = \begin{pmatrix} \frac{E_1^{UD}}{1 - \nu_{12}^{UD} \nu_{21}^{UD}} & \frac{\nu_{12}^{UD} E_2^{UD}}{1 - \nu_{12}^{UD} \nu_{21}^{UD}} & 0 \\ \frac{\nu_{12}^{UD} E_2^{UD}}{1 - \nu_{12}^{UD} \nu_{21}^{UD}} & \frac{E_2^{UD}}{1 - \nu_{12}^{UD} \nu_{21}^{UD}} & 0 \\ 0 & 0 & G_{12}^{UD} \end{pmatrix}. \qquad (55)$$

Table 3 gives the results for both materials. The elements of the honeycomb matrix $Q_H^*$ are significantly lower than the elements of the unidirectional matrix $Q_{UD}^*$.

Now, the results for the total bending rigidity and transverse shear matrix depend on the each ply of the multilayer. We will look at this.

## 4.2. Analytical bending matrix calculation

From Eqs. (34) to (38), the bending rigidity elements of the multilayer can be expressed with the elements of the reduced stiffness matrices of the honeycomb and UD. First, we modify the indexation form of the plies presented in Fig. 3. Since the sandwich is symmetric, we only need half of the multilayer. With this situation, the half-central ply of the honeycomb is

<table><caption>**Table 3**
Elements of the material reduced stiffness matrices of UD and honeycomb.</caption><thead><tr><th>In GPa</th><th>{Q<sup>*</sup><sub>ap</sub>}<sub>UD</sub></th><th>{Q<sup>*</sup><sub>ap</sub>}<sub>H</sub></th></tr></thead><tbody><tr><td>Q<sup>*</sup><sub>11</sub></td><td><b>1.343E+02</b></td><td><b>1.421E-03</b></td></tr><tr><td>Q<sup>*</sup><sub>22</sub></td><td><b>7.738E+00</b></td><td><b>1.421E-03</b></td></tr><tr><td>Q<sup>*</sup><sub>12</sub></td><td><b>2.244E+00</b></td><td><b>1.279E-03</b></td></tr><tr><td>Q<sup>*</sup><sub>66</sub></td><td><b>3.100E+00</b></td><td><b>1.620E-04</b></td></tr><tr><td>Q<sup>*</sup><sub>16</sub></td><td>0</td><td>0</td></tr><tr><td>Q<sup>*</sup><sub>26</sub></td><td>0</td><td>0</td></tr></tbody></table>
Classical multilayer nomination

[UD[$\theta_i, \theta_2, ..., \theta_{l/2}$]$_S$, Honeycomb $\theta_0$]$_S$

Symmetric sandwich with an even number of plies $n$

<figure><img src="image_8.png" /><figcaption>Fig. 6. New indexation form of the plies with $r$.</figcaption></figure>

now noted as row $r=0$ and the others above entire UD plies are noted as rows $r=1$ to $l$ (see Fig. 6). The expression of the bending rigidity elements matrix of the sandwich $D_{\alpha\beta}$ is now written with the elements of the reference reduced stiffness matrices $Q_{\alpha\beta}^r$ from each ply:

$$D_{\alpha\beta} = \frac{2}{3} \times \left( \sum_{r=0}^{\text{external ply } l=(n-1)/2} \sum_{\text{central ply}} d_r Q_{\alpha\beta}^r \right). \qquad (56)$$

The number 2 in Eq. (56) represents the two symmetric parts of the sandwich. The terms $d_r$ are the geometric contribution of each ply with row $r$. They are a function simply of the row $r$, the honeycomb thickness $h_H$ and the composite ply thickness $e_p$. It is understood that the UD plies have rigorously identical thicknesses $e_p=e_{UD}$, so the terms $d_r$ may be written as

$$d_{r=0} = \frac{h_H^3}{8}, \quad d_{r=1 \text{ to } l} = 3 \frac{h_H^2}{4} e_{UD} + 3 \frac{h_H}{2} (2r-1) e_{UD}^2 + (3r^2 - 3r + 1) e_{UD}^3. \qquad (57)$$

Taking into account the rotation $\theta_r$ of each ply and the transformation rule, all terms of Eq. (56) can be detailed as the following condensed form Eq. (58), with simple columns $\{D_{\alpha\beta}\}$, $\{Q_{\alpha\beta}^*\}_H$ and $\{Q_{\alpha\beta}^*\}_{UD}$

$$\{D_{\alpha\beta}\}_{SHC} = 2/3 [\mathbf{C}_1 \cdot \{Q_{\alpha\beta}^*\}_{UD} + d_0 \cdot \mathbf{C}_0 \cdot \{Q_{\alpha\beta}^*\}_{H}], \qquad (58)$$

where

$$\mathbf{C}_1 = \begin{pmatrix} \sum_{r=1}^{l} d_r c_{-\theta_r}^4 & \sum_{r=1}^{l} d_r s_{-\theta_r}^4 & 2 \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & 4 \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 \\ \sum_{r=1}^{l} d_r s_{-\theta_r}^4 & \sum_{r=1}^{l} d_r c_{-\theta_r}^4 & 2 \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & 4 \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 \\ \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & \sum_{r=1}^{l} d_r (c_{-\theta_r}^4 + s_{-\theta_r}^4) & -4 \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 \\ \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & -2 \sum_{r=1}^{l} d_r c_{-\theta_r}^2 s_{-\theta_r}^2 & \sum_{r=1}^{l} d_r (c_{-\theta_r}^2 - s_{-\theta_r}^2)^2 \\ \sum_{r=1}^{l} d_r c_{-\theta_r}^3 s_{-\theta_r} & - \sum_{r=1}^{l} d_r c_{-\theta_r} s_{-\theta_r}^3 & \sum_{r=1}^{l} d_r c_{-\theta_r} s_{-\theta_r} (s_{-\theta_r}^2 - c_{-\theta_r}^2) & 2 \sum_{r=1}^{l} d_r c_{-\theta_r} s_{-\theta_r} (s_{-\theta_r}^2 - c_{-\theta_r}^2) \\ \sum_{r=1}^{l} d_r c_{-\theta_r} s_{-\theta_r}^3 & - \sum_{r=1}^{l} d_r c_{-\theta_r}^3 s_{-\theta_r} & \sum_{r=1}^{l} d_r c_{-\theta_r} s_{-\theta_r} (c_{-\theta_r}^2 - s_{-\theta_r}^2) & 2 \sum_{r=1}^{l} d_r c_{-\theta_r} s_{-\theta_r} (c_{-\theta_r}^2 - s_{-\theta_r}^2) \end{pmatrix}. \qquad (59)$$
Product $d_0 \cdot C_0$ is similar to matrix $C_1$, but is used only for the central half-honeycomb ply. In this case, the index r is limited only to the value 1. Since $d_r$ equals $d_0$ anywhere, and because the angle of honeycomb $\theta_0$ is equal to zero, $C_0$ is simply

$$ C_0 = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. \qquad (60) $$

Compared to {$Q_{\alpha\beta}$)$_H$ and {$Q_{\alpha\beta}$)$_UD$ which have six elements, {$Q_{\alpha\beta}^*$)$_H$ and {$Q_{\alpha\beta}^*$)$_UD$ are reduced to the first four elements. It is understood that $Q_{16}^*$ and $Q_{26}^*$ always tend to zero, which is true for the honeycomb and UD plies (see Table 3). Each element $D_{\alpha\beta}$ of the bending matrix $D^f$ becomes a linear function of the elements of the material reduced stiffness matrices of the honeycomb and UD ($Q_H^*, Q_{UD}^*$).

Note that if all composite plies are replaced (ply per ply) by a unique isotropic material, such as $e=l.e_p$, the previous complex Eq. (58) is equivalent to Eq. (8). Since the mechanical properties of the honeycomb are much less than those of the composite skins, the second term of Eq. (58) can often be neglected, yielding the following expression:

$$ \{D_{\alpha\beta}\}_{SHC}^{\text{simplified}} \approx (2/3)C_1 \cdot \{Q_{\alpha\beta}^*\}_{UD} \qquad (61) $$

This situation is similar to that applied to the bending stiffness of the honeycomb sandwich with isotropic skins. This result is interesting from two points of view: the simple linear function and the global properties of the bending matrix. With the knowledge of the experimental elements $D_{\alpha\beta}$, the first point gives the possibility of coming up to the properties of the UD (application at Appendix A). The second point is to control the final orthotropic properties between the elements of the bending matrix, which is going to depend mainly on the combination of laminate angles.

### 4.3. Analytical transverse shear matrix calculation

With the laminate shown in Fig. 3, one starts out again from Eqs. (39) to (43). Berthelot [12] has developed formulations for $C_{11}^{\text{hom}}$, $C_{22}^{\text{hom}}$ and $C_{12}^{\text{hom}}$ for a laminate. It is assumed that all plies of UD have the same thickness or $e_k=e_{UD}$, for whatever k. The material axes (1,2,3) of the UD at 0° and (L,T,T) of the honeycomb are parallel to the reference axes (x,y,z). Thus, we have

$$ C_{11}^{\text{hom}} = k_{11} \cdot \frac{\left[ e_{UD} \cdot \sum_{k=1 \text{ to } n}^{\text{(except the honeycomb)}} (G_{23}^{\text{UD}} \sin^2 \theta_k + G_{13}^{\text{UD}} \cos^2 \theta_k) \right] + h_H G_{xz}^H}{(h_H + 2e)}, \qquad (62) $$

$$ C_{22}^{\text{hom}} = k_{22} \cdot \frac{\left[ e_{UD} \cdot \sum_{k=1 \text{ to } n}^{\text{(except the honeycomb)}} (G_{23}^{\text{UD}} \cos^2 \theta_k + G_{13}^{\text{UD}} \sin^2 \theta_k) \right] + h_H G_{yz}^H}{(h_H + 2e)}, \qquad (63) $$

$$ C_{12}^{\text{hom}} = k_{12} \cdot \frac{\left[ e_{UD} \cdot \sum_{k=1 \text{ to } n}^{\text{(except the honeycomb)}} (G_{13}^{\text{UD}} - G_{23}^{\text{UD}}) \cos \theta_k \sin \theta_k \right] + 0}{(h_H + 2e)}, \qquad (64) $$

For our application, one symmetric face laminate is composed $(0^\circ, +45^\circ, +90^\circ, -45^\circ)_s$. Also, for a half-face, the trigonometric functions give $\Sigma_{\theta=0.45,90,-45} \cos^2 \theta_k = 2$, $\Sigma_{\theta=0.45,90,-45} \sin^2 \theta_k = 2$, $\Sigma_{\theta=0.45,90,-45} \cos \theta_k \sin \theta_k = 0$.

The two faces represent four half-faces and the elements of the $C_{11}^{\text{hom}}$ matrix are calculated in such situations as follows:

$$ C_{11}^{\text{hom}} = k_{11} \cdot \frac{8e_{UD}G_{13}^{\text{UD}} + 8e_{UD}G_{23}^{\text{UD}} + h_H G_{xz}^H}{(h_H + 2e)} = k_{11} \cdot \frac{e(G_{13}^{\text{UD}} + G_{23}^{\text{UD}}) + h_H G_{xz}^H}{(h_H + 2e)}, \qquad (65) $$

$$ C_{22}^{\text{hom}} = k_{22} \cdot \frac{8e_{UD}G_{13}^{\text{UD}} + 8e_{UD}G_{23}^{\text{UD}} + h_H G_{yz}^H}{(h_H + 2e)} = k_{22} \cdot \frac{e(G_{13}^{\text{UD}} + G_{23}^{\text{UD}}) + h_H G_{yz}^H}{(h_H + 2e)}, \qquad (66) $$

$$ C_{12}^{\text{hom}} = 0. \qquad (67) $$

If the composite UD becomes an isotropic material, it is as $C_{13}^{\text{UD}} = G_{23}^{\text{UD}} = G$. We again obtain the formulas of Eq. (11), if we associate $k_{11}$ and $k_{22}$ of Berthelot [12] with $1/k_{xz}$ and $1/k_{yz}$ of Gay [11]. It is the result for a symmetric sandwich with isotropic faces.

The analytical $k_{xz}$ and $k_{yz}$ developments are supplied by Gay for a three-layer case, but with a complex sandwich, obtaining the analytical formulas for $k_{11}$ and $k_{22}$ is not obvious. To avoid this difficulty, we must know $C_{11}^{\text{hom}}$ and $C_{22}^{\text{hom}}$ directly. This can be achieved with the help of MSC/Nastran. With a CQUAD4 card associated to a PCOMP card and asking
(PARAM, PRTPCOMP, YES), we obtain an equivalent PSHELL card in the PUNCH file. This PSHELL refers to three news MAT2 cards. These three MAT2 give the elements of the homogenized membrane A<sup>hom</sup>, flexural D<sup>hom</sup> and transverse shear C<sup>hom</sup> matrices:

$$A^{\text{hom}} = \frac{1}{h_s} \cdot A^M, \qquad (68)$$

$$D^{\text{hom}} = \frac{12}{h_s^3} \cdot D^F, \qquad (69)$$

$$C^{\text{hom}} = \frac{1}{h_s} \cdot C^S, \qquad (70)$$

where, for the elements of the matrices A<sup>M</sup>, D<sup>F</sup> and C<sup>S</sup>:

$$A_{\alpha\beta} = \int_{-(h_s/2)}^{h_s/2} Q_{\alpha\beta} dz = \sum_{k=1}^{n} \int_{z_{k-1}}^{z_k} Q_{\alpha\beta}^k dz = \sum_{k=1}^{n} (z_k - z_{k-1}) Q_{\alpha\beta}^k, \qquad (71)$$

$$D_{\alpha\beta} = \int_{-h_s/2}^{+h_s/2} z^2 Q_{\alpha\beta} dz = \sum_{k=1}^{n} \int_{z_{k-1}}^{z_k} z^2 Q_{\alpha\beta}^k dz = \frac{1}{3} \sum_{k=1}^{n} (z_k^3 - z_{k-1}^3) Q_{\alpha\beta}^k$$

previously defined Eq. (34),

$$C_{11} = k_{11} \sum_{k=1}^{n} e_k G_{xz}^k, \qquad (72)$$

$$C_{22} = k_{11} \sum_{k=1}^{n} e_k G_{yz}^k, \qquad (73)$$

C<sub>11</sub> and C<sub>22</sub> given by Nastran are assumed to be condensed format numerators of C<sub>11</sub><sup>hom</sup> Eq. (41) and C<sub>22</sub><sup>hom</sup> Eq. (42). Finally, S for our application is estimated in the same way as Eq. (39), in which C<sub>11</sub><sup>hom</sup> and C<sub>22</sub><sup>hom</sup> are given by the MSC/Nastran software:

$$S_{\text{Nastran}} = h_s (C_{11\text{Nast}}^{\text{hom}} \cdot C_{22\text{Nast}}^{\text{hom}})^{1/2} = (C_{11\text{Nast}} \cdot C_{22\text{Nast}})^{1/2}. \qquad (74)$$

C<sub>11</sub><sup>hom</sup> and C<sub>22</sub><sup>hom</sup> from the third MAT2 are no more no less than the homogenized transverse shear moduli in the reference planes xz and yz for this sandwich with composites faces, as G<sub>xzeq</sub> and G<sub>xyeq</sub> for the symmetric honeycomb sandwich with isotropic faces into Eq. (13) S = h<sub>s</sub>(G<sub>xzeq</sub> · G<sub>xyeq</sub>)<sup>1/2</sup>.

The MSC/Nastran software gives us for our own application:

$$C_{11}^{\text{hom}} = 50788708 \, \text{N/m}^2, \quad C_{22}^{\text{hom}} = 33350626 \, \text{N/m}^2, \quad C_{12}^{\text{hom}} = 214 \, \text{N/m}^2, \quad C_{13,23,33}^{\text{hom}} = 0, \\ \text{whether with } h_s = 13.910^{-3} \, \text{m}, \quad C_{11} = 705963 \, \text{N/m}, \quad C_{22} = 463573 \, \text{N/m} \text{ and } S = 572071 \, \text{N/m}.$$

<figure><img src="image_9.png" /><figcaption>Fig. 7. Comparison between the real and estimated geometric transverse shear moduli of the honeycomb as a function of h<sub>H</sub>.</figcaption></figure>
It is interesting to see the influence of the thickness of the honeycomb for this type of sandwich, and compare it to a sandwich with isotropic skins. In Fig. 7, the transverse shear stiffness of the sandwich is plotted as a function of the thickness $h_N$ of the honeycomb, while the face thickness is maintained constant at 1 mm. S is computed with Nastran from e (74) $S_{\text{Numerical}} = h_S (C_{11Nast}^{\text{hom}} C_{22Nast}^{\text{hom}})^{1/2}$ and with analytical equation (13) $S_{\text{Analytical}} = h_S (G_{xzeq}^{\text{Allen}} G_{yzeq}^{\text{Allen}})^{1/2}$ in which Allen transverse shear moduli equation (14) are used. $G_{xzeq}^{\text{Allen}}$ and $G_{yzeq}^{\text{Allen}}$ are $G_{xzeq}^{\text{Simpl}}$ and $G_{yzeq}^{\text{Simpl}}$ of the symmetric honeycomb with isotropic faces with particular conditions: $h_H \gg 2e$, $G_{\text{face}} \gg G_{xz}^H$ and $G_{\text{face}} \gg G_{yz}^H$. In the graph (Fig. 7) and for the sandwich with composites faces, there is no difference between the two methods for $h_H$ greater than 1.25 mm. It means in this zone, $C_{11Nast}^{\text{hom}} = G_{xzeq}^{\text{Allen}}$ and $C_{22Nast}^{\text{hom}} = G_{yzeq}^{\text{Allen}}$. We are in the same situation as that of the symmetric honeycomb with isotropic faces. With a small face thickness, S is once again only a function of e, $h_H$, $G_{xz}^H$ and $G_{yz}^H$: $S_{\text{Analytical}} = h_H (1+(e/h_H))^2 G_g^H$ with $G_g^H = (G_{xz}^H G_{yz}^H)^{1/2}$.

The skins have no influence except the thickness. Conversely, we can estimate the geometric transverse shear module of the honeycomb from S in the following manner:

$$G_g^H \text{ estimated} = (G_{xz}^H G_{yz}^H)^{1/2} \text{ estimated} = \frac{h_H}{(h_H + e)^2} S. \quad (75)$$

This is confirmed by the curve of Eq. (75) in Fig. 7, in which $S_{\text{Analytical}}^{\text{Allen}}$ is substituted by $S_{\text{Numerical}}^{\text{Nastran}}$. The estimated value of the geometric transverse shear modulus of the honeycomb converges to the real value.

### 4.3.1. Numerical application

With $e=1$ mm, $h_H=11.9$ mm and $S_{\text{Nastran}}=572,071$ N/m, the estimated geometric transverse shear modulus of the nomex honeycomb is $4.09E+07$ N/m$^2$. The real value is $4.08E+07$ N/m$^2$.

Finally, though we are able to estimate S through the identification of the wavenumber or modal density or mechanical impedance for example, the geometric transverse shear modulus of the honeycomb can also be estimated from experiments.

### 4.4. Analytical mass per unit area

The Nomex® honeycomb of 12 mm thickness $h_H$ has a density $\rho_H$ of 53 kg/m$^3$. The two composite faces have a thickness $e$ of 1 mm and a density $\rho_{UD}$ equal to 1550 kg/m$^3$. The theoretical mass per unit area is 3.73 kg/m$^2$. The common problem is to take into account the glue placed at the honeycomb skin interfaces and sometimes densification with light foam included at the contour of the sandwich. We can easily estimate this by weighing the panels. This additional mass finally gives a real mass per unit area of 5.83 kg/m$^2$ for the first panel A of dimensions 1.96(L) x 1.36(T) m$^2$ and 4.68 kg/m$^2$ for the second panel B of dimensions 1(L) x 1(T) m$^2$. The two panels seem to be identical multilayer structures, but the mass per unit area dispersion comes from the local distribution of the glue and the densification. The modeling of this aspect is not perfect. During numerical computation and in order to take this aspect into account, we artificially increased the density of the honeycomb to find the same surface density. The new total density (honeycomb+glue+densification) has been set respectively to 230 kg/m$^3$ and 133 kg/m$^3$.

### 4.5. Analytical transition and critical frequencies

Table 4 summarizes the characteristics of the two panels tested. The transition and critical frequencies are calculated with $D=(D_{11},D_{22})^{1/2}$. Fig. 8 locates these frequencies on the analytical curve of the wavenumber. We use the simplified equation (20) for the bending wavenumber of the panels and Eq. (30) for the acoustic wavenumber.

#### 4.5.1. Analytical transition frequency

The main properties of our own symmetric honeycomb panels with composite faces are already defined: $D_{11}=4621$ N m, $D_{22}=4611$ N m, $D_{12}=1490$ N m, $D_{66}=1561$ N m, $S=572,071$ N/m, $m_{SA}=5.83$ kg/m$^2$ and $m_{SB}=4.68$ kg/m$^2$. The transition frequency can be estimated with Eq. (29). Since $D_{11} \approx D_{22} \approx H \approx (D_{11},D_{22})^{1/2}$, the complex sandwich has a bending stiffness

<table><caption>**Table 4**
Analytical transition and critical frequencies.</caption><thead><tr><th>Panel</th><th>Dimensions (m) × (m)</th><th>S (N/m)</th><th>m<sub>S</sub> (kg/m<sup>2</sup>)</th><th>D<sub>11</sub> (N m)</th><th>D<sub>22</sub> (N m)</th><th>D (D<sub>11</sub>,D<sub>22</sub>)<sup>1/2</sup> (N m)</th><th>H (D<sub>12</sub>+2D<sub>66</sub>) (N m)</th><th>f<sub>trans</sub>(m<sub>S</sub>,D,S) transition frequency (Hz)</th><th>Criterion to compute f<sub>cr</sub> (1 - m<sub>S</sub>c<sub>air</sub><sup>2</sup>/S) ≥ 0</th><th>f<sub>cr</sub>(m<sub>S</sub>,D,S) critical frequency (Hz)</th></tr></thead><tbody><tr><td>A</td><td>1.96(L) × 1.36(T)</td><td>572,071</td><td>5.83</td><td>4621</td><td>4611</td><td>4616</td><td>4612</td><td>278</td><td>-0.178</td><td>Indetermination</td></tr><tr><td>B</td><td>1(L) × 1(T)</td><td>572,071</td><td>4.68</td><td>4621</td><td>4611</td><td>4616</td><td>4612</td><td>310</td><td>0.053</td><td>2534</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Analytical transition and critical frequency positions of panels A and B in the frequency domain.</figcaption></figure>

behavior that is analogous to that of an isotropic panel. *D* is taken as the geometric stiffness of D₁₁ and D₂₂ either 4616 N m. The transition frequencies of Panels A and B are respectively 278 Hz and 310 Hz. They are very low.

### 4.5.2. Analytical critical frequency

The critical frequency is estimated with Eq. (31). With $c_{air}=340$ m/s, the critical frequency of Panel A is undefined. The critical frequency of Panel B is 2534 Hz. For both panels, which have the same material constitution, *D* and *S* are invariant. The only difference is $m_s$. The mass per unit area of the first panel $m_{SA}$ yields a negative value in the denominator of the critical frequency. From a mathematical point of view, the critical frequency cannot be known. The bending wavenumber modulus is higher than the acoustic wavenumber. It is as if the intersection of the bending wavenumber modulus with $k_{ac}$ was rejected at a higher frequency. The difference of modification between the two panels, A and B, is the quantity of glue and densification around the contour of the sandwich. The problem is for the experimental and theoretical mass per unit area to agree. Given the difficulty in obtaining the true theoretical value of $m_s$, the critical frequency of the first panel A is approximately estimated between 2500 Hz and a greater frequency. Another possible explanation is the model used. The analytical bending wavenumber formula depends on the simplifications made during the presentation of the differential equation of displacement. The original model with a 6th order differential equation is reduced to the fourth order. If the bending wavenumber was described by a 6th order differential equation, the model would be more accurate and the module amplitude will go down in the high frequencies. Consequently the exact solution will, in the frequency domain, always result in an interception with the wavenumber in air, giving a critical frequency. We can see this in the following analysis, done with experimental and numerical developments.

# 5. Experimental, analytical and numerical results

## 5.1. Numerical model description

In order to evaluate all of the physical quantities (wavenumber, frequencies, modal density and mechanical impedance) of the complex sandwich with the only three parameters ($m_s$, *D*, *S*), numerical models are developed with MSC/Nastran. Although $h_H$ equals 14 mm, the sandwich panel is still considered thin with transverse shear. The criterion for isotropic plate $L_{min}(a,b)/h_H$ is well above 20. The field of study is rather between the Kirchhoff-Love and Reisner-Mindlin theories. The sandwich is modeled as one multilayer. The honeycomb core and composite plies properties are fully included through a PCOMP card. The surfaces of the rectangular and square plates are meshed with the CQUAD4 finite element. The dimensions of the surface mesh size are 1 cm by 1 cm. The model has a single layer of CQUAD elements. The total thickness is represented by the multilayer PCOMP card. The pair of cards CQUAD and PCOMP automatically generates a PSHELL card. At the end, this new card contains all properties: membrane, bending and out-of-plane shear. This model is called the "2D model".

With MSC/Nastran, in terms of the quantities sought, various computations are requested: MAT2 of PSHELL properties in the Punch file, the eigenfrequencies and the deformation of the modal shapes in the F06 file. For the wavenumber and mechanical impedance characterizations, we include excitation points. In this case, frequency response functions (FRF) for some points are required. The total loss factor of the PCOMP is taken as 0.014 over the entire frequency band analyzed. Different boundary conditions around the plate are also used: free, simply supported and clamped.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Description of the measurement area.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Rotations $\psi$ and dimensions $r_a$ of the antenna.</figcaption></figure>

For the free boundary conditions, there is no particular condition around the plate. In the case of the simply supported boundary conditions we have

- along the two edges a: the displacements $u_x=v_y=w_z=0$ and the rotations $R_y=R_z=0$,

- along the two edges b: the displacements $u_x=v_y=w_z=0$ and the rotations $R_x=R_z=0$.

The clamped boundary conditions are along all the edges: the displacements $u_x=v_y=w_z=0$ and the rotations $R_x=R_y=R_z=0$.

## 5.2. Wavenumber modulus under free and clamped boundary conditions

Experimental, analytical and numerical wavenumber moduli are compared in Fig. 12 for the first panel A, under free boundary conditions. The analytical wavenumber modulus could be calculated from Eq. (54), but since the orthotropic coefficient $H/D$ is close to unity, it is directly done from the simplified Eq. (20) $|k(\omega)| = ((1/2S)(m_s\omega^2 + \omega(m_s^2\omega^2 + (4m_sS^2/D))^{1/2}))^{1/2}$. Experimental and numerical wavenumber moduli are strictly established from the same following formulation (76) for a thin plate. This was proposed by Carniel and Gavric [21] for the measurement of the structural intensity and the software developed to compute it by a Finite Element Method:

$$|k(x_m, y_m, f)| = \sqrt{\frac{4S_{11,xm,ym} - \text{Re}(S_{12,xm,ym} + S_{13,xm,ym} + S_{14,xm,ym} + S_{15,xm,ym})}{r_a^2 S_{11,xm,ym}}} \quad (76)$$

The panel is excited by a mechanical force at a driving point. We define the measurement area (see Fig. 9) as a small part of the panel whose center is the point $P_1(x_m, y_m)$. From this central point, four points ($P_2-P_5$) are chosen on a same circle whose radius is $r_a$. The point distribution forms an orthogonal cross. $S_{11}$ is the displacement autopower spectrum at central
<figure><img src="image_13.png" /><figcaption>Fig. 11. Locations of the driving point and zones of Panel A, for $k(f)$ and $Z_\infty(f)$ analysis.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Analytical, numerical and experimental wavenumber modulus $k(f)$ of panel A, under free boundary conditions, (Zone 3, $r_a=2$ cm).</figcaption></figure>

point $P_1, S_{12}, S_{13}, S_{14}$ and $S_{15}$ are the displacement cross-power spectra of points $P_2, P_3, P_4$ and $P_5$ related to point $P_1$. This antenna is able to turn around the central point to obtain an angle $\psi$ (Fig. 10). The possibilities of the antenna rotation give measurements as a function of this angle, in order to analyze whether the material is isotropic or orthotropic. Radius $r$ is adapted to the bandwidth frequency as a function of the wavelength.

Vibration displacement measurements on Panel A were carried out with a laser velocimeter device. Several zones of the panel have been used for the wavenumber characterization (Fig. 11).

For Panel A in Zone 3 with an antenna radius of 2 cm, Fig. 12 shows a good agreement between experimental, analytical and numerical results. The transition frequency is located at 278 Hz. The experimental results are rather independent for the three angles $\psi$ tested. The analytical and numerical results are above that of the acoustical wavenumber, while the experiment results tend to cross that of acoustics at high frequencies. The analytical and numerical curves follow a straight line having a common slope of $2\pi(m_s/S)^{1/2}$. The acoustical wavenumber slope is $2\pi/c_{\text{air}}$. The critical frequency is not obvious to locate. For analytical and numerical results, it tends to infinity and for experimental results it seems to be around 3500 Hz.

Extensive research shows that if we change, for example, the area or the boundary conditions, we do not change the numerical modulus of the wavenumber. In Fig. 13, the panel is at first under free and then clamped boundary conditions. The slope of the wavenumber modulus is the same over the transition frequency. The wavenumber is fully set to the intrinsic material properties of the sandwich. In the high frequency domain, it is mainly caused by mass per unit area $m_s$ and transverse shear $S$.

On the other side, experiments and numerical modeling of the wavenumber depend on the choice of the value attributed to the antenna radius $r_a$. With a larger radius of 5 cm and still the same number of modes requested for numerical computation (1800), the levels of the wavenumber modulus (experiments and numerical) decrease after 3000 Hz. This is due to a known issue with the ratio between the antenna radius and the analyzed wavelength. With an antenna radius of 2 cm and a sufficiently high mode number, the analytical theory and numerical modeling of the wavenumber are rather in phase. If the mass per unit area is reduced without changing the transverse shear $S$, we will gradually reach the wavenumber of Panel B, as in Fig. 8.

After the analysis of the wavenumber, the effect of the transverse shear stiffness $S$ is now observed on the modal frequencies.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Numerical wavenumber modulus *k*(f) in the four zones of Panel A ($r_a=2$ cm), under free then clamped boundary conditions.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Analytical and numerical modal eigenfrequencies of Panel A, under simply supported boundary conditions.</figcaption></figure>

### 5.3. Modal frequencies under the simply supported boundary conditions

The wavenumber in the sandwich is clearly analyzed and controlled; we will look at the impact on the modal frequencies of the larger panel. Panel A is only analyzed under simply supported boundary conditions. The analytical modal frequencies of Panel A are computed with Eq. (50) and then compared to the first 50 numerical frequencies obtained from MSC/Nastran, in Fig. 14. The rows at the bottom of the graph are sorted in ascending order of the associated analytical frequencies. The rows of the reclassed numerical modes (*i*,*j*) are according to the analytical rows (*i*,*j*).

The numerical frequencies are in good agreement with the analytical frequencies. The difference between the 100 first frequencies gradually increases up to 9 percent. All physical quantities ($m_s$, $D_{αβ}$, $S$) are well taken into account. From these frequencies, it is easy to analyze the modal density of Panel A in the next section.

### 5.4. Modal density under the simply supported boundary conditions

The modal density of Panel A is computed from the previous analytical frequencies and compared to the analytical modal density from Eq. (26), in Fig. 15. A sliding window along the frequency axis is used. We calculate for each central frequency of the window, the average spread between plus or minus five neighboring frequencies. Then, the average spread is inversed to obtain the mean density for the central frequency. The stack of one face [$0°, 45°, 90°, -45°]$ is quasi-isotropic composite.
<figure><img src="image_17.png" /><figcaption>Fig. 15. Analytical modal density of Panel A, under simply supported boundary conditions.</figcaption></figure>

The modal density of two faces [[0°,45°, 90°, -45°]s]₂ (or the sandwich without honeycomb) is also plotted. The initial analytical frequencies are computed with the Nastran properties from the PUNCH output file.

Properties of the stack with two skins [[0°, 45°, 90°, -45°]s]₂

$$D_{11} = 42.82 \text{ N m}, \quad D_{22} = 32.94 \text{ Nm}, \quad D_{12} = 10.94 \text{ N m}, \quad D_{66} = 11.51 \text{ N m}, \quad H = 33.97 \text{ N m},$$

$$\alpha = H/D_{11} = 0.8, \quad S = 4630126 \text{ N/m} \quad \text{and} \quad m_{S2e} = 3.1 \text{ kg/m}^2.$$

The transition frequencies of the two skins together and Panel A are respectively estimated at 31,979 Hz and 278 Hz. Since the transition frequency of the faces gathered is very high, the modal density is like that of a thin isotropic plate. The analytical points of the modal density are around a constant mean value, and the modal density is rather independent of the frequency. The introduction of the honeycomb between the two faces increases the mass per unit $m_s$ because the glue, especially increases the stiffness matrix $D_{\alpha\beta}$ and decreases the transverse shear $S$ for the sandwich. The analytical modal density of Panel A at a low frequency is lower than that of the two faces together due to the increase in stiffness. Thus, above the transition frequency, the modal density begins to be more sensitive to transverse shear. The modal density constant at a low frequency, now increases with the frequency. The bending modal frequencies are increasingly tight. Such a behavior is totally comparable to a honeycomb sandwich with isotropic skins.

## 5.5. Mechanical impedance under the free boundary conditions

After having analyzed the wavenumber modulus and modal density, it is interesting to observe the transition frequency on the mechanical impedance for such a honeycomb sandwich. The driving point selected on Panel A is marked in Fig. 11. The inverse of the mechanical impedance at the excitation point is computed from the finite element model with MSC/Nastran. The real part is conserved and then compared to the complete analytical formulation equation (28). The numerical and analytical results are plotted in Fig. 16.

The total loss factor introduced in the PCOMP card is 0.014 for the entire frequency band analyzed. This stems from experimental results. Once again, the transition frequency separates the band between two zones. At low frequencies, the sandwich has an isotropic plate behavior. Above the transition frequency, the average numerical curve is perfectly on the analytical curve. The analysis is the same for Panel B. The driving point selected on Panel B is marked in Fig. 17. The inverse of the mechanical impedance at the excitation point is plotted in Fig. 18. The stiffness matrix elements $D_{\alpha\beta}$ and $S$ are identical for both panels. Only the mass per unit area $m_s$ makes the difference for the global material of the multilayer. For this reason and from Table 4, the transition frequency is estimated at 310 Hz. Panels A and B have different dimensions. The frequency distribution of the modes changes but not the slope of the inverse of the mechanical impedance at high frequencies. At low frequencies, if we want the numerical result to agree perfectly with the analytical result, it is necessary to obtain more driving points. The SEA formulations need to have a great number of modes and this is not the case at low frequencies. Another faster method is to integrate the amplitude response over a large bandwidth frequency or to greatly increase the damping.
<figure><img src="image_18.png" /><figcaption>Fig. 16. Analytical and numerical inverse of the mechanical impedance, at driving point of Panel A, under free boundary conditions.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Location of the driving point of Panel B for $Z_\infty(f)$ analysis.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Analytical and numerical inverse of the mechanical impedance, at the driving point of the square panel B, under free boundary conditions.</figcaption></figure>
## 6. Presentation of a formulation of modal frequencies under the free boundary conditions

With all parameters being controlled simultaneously with analytical and numerical operations, it is possible to identify some honeycomb and composite properties from dynamic measurement and through a special process. At the end of paragraphs 4.2 and 4.3 we show how we can trace back by simulation, on the one hand to the physical properties of the unidirectional composite ($E_1^{UD}$, $E_2^{UD}$, $G_{12}^{UD}$ and $\nu_{12}^{UD}$), and on the other hand to the geometric transverse shear modulus of the honeycomb $G_g^{HD}$. In order to return to these, it is necessary to first determine $D_{11}$, $D_{22}$, $D_{12}$, $D_{66}$ and S. From the experimental point of view, S can be estimated by simply measuring the mechanical impedance, or the modal density, or the modulus of the wavenumber, which is more complicated. The coefficients $D_{\alpha\beta j}$ of the bending matrix need the modal frequencies of an experimental modal analysis. The coefficients are in theory independent from the boundary conditions. It is not easy to obtain real simply supported or clamped boundary conditions during experiments. Since the frequency values often deviate from the theoretical frequencies, the free boundary conditions are preferred. Under the conditions of this case, we must develop a formula that takes into account the stiffness matrix and the out-of-plane shear.

### 6.1. Methodology and development of the formula for free boundary conditions

Looking again at Eq. (50) for a simply supported plate, we observe that the first term between the large brackets represents the modal frequencies without transverse shear. The second term transforms each frequency without transverse shear into a new frequency taking into account the shear stiffness S. This is the effect of the presence of the transverse shear. We can easily imagine that this effect also exists for a free or clamped plate. Of course, the correction factor of the out-of-plane shear is not necessarily the same. We call it $C_{\varphi, i,j}^S$, where $\varphi$ represents the boundary conditions (free, or simply supported, or clamped) and $(i,j)$ are the rows of the mode. The frequencies without transverse shear are well known for an isotropic and orthotropic plate as a function of the boundary conditions. Warburton [14] was the first to give the formulations for an isotropic plate using $(m,n)$ rows in terms of a stationary sine wave and later Blevins [7] studied them for an orthotropic plate using $(i,j)$ rows in terms of a half-stationary sine wave. The generalized modal frequency formula can be written in a general form with the indices $(m,n)$ or the indices $(i,j)$, as follows:

with Warburton annotations:

$$f_{\varphi,m,n}^{S} = \frac{\pi}{2a^2} \left( \frac{\langle \lambda_{\varphi,m,n}^{2\text{Warburton}} ; D \rangle}{m_s} \right)^{1/2} \cdot C_{\varphi,m,n}^{S}, \quad (77)$$

and

$$\left\{ \begin{aligned} \langle \lambda_{\varphi,m,n}^{2\text{Warburton}} ; D \rangle &= G_{x_{\varphi,m}}^4 + G_{y_{\varphi,n}}^4 \frac{a^4}{b^4} + 2 \frac{a^2}{b^2} (\nu H_{x_{\varphi,m}} H_{y_{\varphi,n}} + (1-\nu)j_{x_{\varphi,m}} j_{y_{\varphi,n}}) \\ \langle \lambda_{\varphi,m,n}^{2\text{Warburton}} ; D \rangle &= a^4 \left( \frac{G_{x_{\varphi,m}}^4}{a^4} + \frac{G_{y_{\varphi,n}}^4}{b^4} + \frac{2j_{x_{\varphi,m}} j_{y_{\varphi,n}} + 2\nu (H_{x_{\varphi,m}} H_{y_{\varphi,n}} - j_{x_{\varphi,m}} j_{y_{\varphi,n}})}{a^2 b^2} \right) D, \end{aligned} \right. \quad (78)$$

$G_{x_{\varphi,m}}$, $G_{y_{\varphi,n}}$, $H_{x_{\varphi,m}}$, $H_{y_{\varphi,n}}$, $j_{x_{\varphi,m}}$, $j_{y_{\varphi,n}}$ are the Warburton analytical coefficients, as a function of $\varphi$, $m$ and $n$, with Blevins annotations:

$$f_{\varphi,i,j}^{S} = \frac{1}{2\pi a^2} \left( \frac{\langle \lambda_{\varphi,i,j}^{4\text{Blevins}} ; D_{\alpha\beta} \rangle}{m_s} \right)^{1/2} \cdot C_{\varphi,i,j}^{S}, \quad (79)$$

and

$$\langle \lambda_{\varphi,i,j}^{4\text{Blevins}} ; D_{\alpha\beta} \rangle = \pi^4 a^4 \frac{G_{1_{\varphi,j}}^4}{a^4 D_{11}} + \frac{G_{2_{\varphi,j}}^4}{b^4} D_{22} + \frac{2}{a^2 b^2} H_{1_{\varphi,j}} H_{2_{\varphi,j}} (D_{12} + 2D_{66}) \\ + \frac{4}{a^2 b^2} D_{66} (j_{1_{\varphi,j}} j_{2_{\varphi,j}} - H_{1_{\varphi,j}} H_{2_{\varphi,j}}) \quad (80)$$

$G_{1_{\varphi,j}}$, $G_{2_{\varphi,j}}$, $H_{1_{\varphi,j}}$, $H_{2_{\varphi,j}}$, $j_{1_{\varphi,j}}$, $j_{2_{\varphi,j}}$ are the Blevins analytical coefficients, as a function of $\varphi$, $i$ and $j$.

Ranks $(m,n)$ are connected to rows $(i,j)$ as a function of the boundary conditions. For the simply supported conditions we have $i=m-1$ and $j=n-1$, while for free conditions we have $i=m+1$ and $j=n+1$. The Blevins analytical coefficients and the Warburton analytical coefficients give one-to-one the same values. When the Blevins orthotropic plate becomes isotropic,
this means that $D_{11} = D_{22} = H = D$, $D_{12} = \nu_{21}D_{11} = \nu_{21}D$ and $D_{66} = (H - \nu_{21}D_{11})/2 = (1 - \nu_{21})D/2$ and we observe

$$ \begin{aligned} \langle \lambda_{\varphi, i,j}^{4\text{Blevins}}; D_{\alpha\beta} \rangle &= \lambda_{\varphi, i,j}^{4\text{Blevins}} \cdot D \quad \text{with} \quad \lambda_{\varphi, i,j}^{4\text{Blevins}} \\ &= \pi^4 \left( G_{1,\varphi,i}^4 + G_{2,\varphi,j}^4 \frac{a^4}{b^4} + 2 \frac{a^2}{b^2} (\nu_{21}H_{1,\varphi,i}H_{2,\varphi,j} + (1-\nu_{21})j_{1,\varphi,i}j_{2,\varphi,j}) \right). \end{aligned} $$

Substituting $\nu_{21}$ for $\nu$, as the pair $(i,j)$ corresponds to the pair $(m,n)$, that is, $\lambda_{\varphi,m,n}^{4\text{Blevins}} = \lambda_{\varphi,i,j}^{2\text{Blevins}}/\pi^2$. There is no difference between $f_{\varphi,i,j}^S$ and $f_{\varphi,m,n}^S$, since $\varphi$ is the same boundary conditions.

For $\varphi$ in the simply supported plate (SSP), $C_{\varphi,i,j}^S$ and $C_{\varphi,m,n}^S$ are

$$ C_{SSP,i,j}^{S} = 1 + \frac{\pi^2}{S} \left( D_{11} \left( \frac{i}{a} \right)^2 + D_{22} \left( \frac{j}{b} \right)^2 \right)^{-1/2}, \quad C_{SSP,m,n}^{S} = \left( 1 + \frac{\pi^2}{S} \left( D_{11} \left( \frac{m-1}{a} \right)^2 + D_{22} \left( \frac{n-1}{b} \right)^2 \right) \right)^{-1/2}. \quad (81) $$

For $\varphi$ in the free plate (FP) and like in the above formula, $C_{\varphi,i,j}^N$ and $C_{\varphi,m,n}^N$ have been tested and the best correlation between analytical and numerical frequencies has been found with

$$ C_{FP,i,j}^{S} = \left( 1 + \frac{\pi^2}{S} \left( D_{11} \left( \frac{i-1}{a} \right)^2 + D_{22} \left( \frac{j-1}{b} \right)^2 \right) \right)^{-1/2}, \quad C_{FP,m,n}^{S} = \left( 1 + \frac{\pi^2}{S} \left( D_{11} \left( \frac{m}{a} \right)^2 + D_{22} \left( \frac{n}{b} \right)^2 \right) \right)^{-1/2}. \quad (82) $$

In order to simplify the presentation and because it is easier to count the number of nodes on mode shape under free conditions, henceforth $(m,n)$ rows are preferably used to $(i,j)$. We introduce $(m,n)$ in Eq. (79), substituting the term $\langle \lambda_{\varphi,i,j}^{4\text{Blevins}}; D_{\alpha\beta} \rangle$ by a new $\langle \lambda_{\varphi,m,n}^{2\text{Warburton}}; D_{\alpha\beta} \rangle$ and we replace $C_{\varphi,i,j}^S$ by $C_{\varphi,m,n}^S$. With the Warburton indices and coefficients, the modal frequencies defined by Blevins become

$$ f_{\varphi,m,n}^{S} = \frac{\pi}{2a^2} \left( \frac{\langle \lambda_{\varphi,m,n}^{2\text{Warburton}}; D_{\alpha\beta} \rangle}{m_s} \right)^{1/2} \cdot C_{\varphi,m,n}^{S}, \quad (83) $$

$$ \langle \lambda_{\varphi,m,n}^{2\text{Warburton}}; D_{\alpha\beta} \rangle = a^4 \left( \frac{G_{x_{\varphi,m}}^4}{a^4} D_{11} + \frac{G_{y_{\varphi,n}}^4}{b^4} D_{22} + \frac{2}{a^2 b^2} H_{x_{\varphi,m}} H_{y_{\varphi,n}} (D_{12} + 2D_{66}) \right. \\ \left. + \frac{4}{a^2 b^2} D_{66} (j_{x_{\varphi,m}} j_{y_{\varphi,n}} - H_{x_{\varphi,m}} H_{y_{\varphi,n}}) \right). \quad (84) $$

In Fig. 19, the modal frequencies of the rectangular plate A under free conditions are plotted using Warburton and Blevins analytical formula, first without and then with transverse shear. The main parameters are for the case of isotropic Plate A $D \approx D_{11} \approx D_{22} \approx H \approx (D_{11}D_{22})^{1/2} = 4616$ N m, $m_{SA}=5.83$ kg/m², $\nu=0.32$, for the case of orthotropic Plate A ($D_{11},D_{22},D_{12},D_{66}, H=(4621,4611,1490,1561,4612)$ N m, $\nu_{21}=\nu_{12}=0.32$, and $S=572,071$ N/m. There is no difference between Warburton and Blevins computations. The numerical result with transverse shear $S$ coming from Nastran software is also included. The analytical and numerical results are in good agreement.

<figure><img src="image_21.png" /><figcaption>Fig. 19. Analytical and numerical modal eigenfrequencies of Panel A, under free boundary conditions.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 20. Analytical and numerical modal eigenfrequencies of the square panel B, under free boundary conditions.</figcaption></figure>

In Fig. 20, the modal frequencies of the square plate B under free conditions are also plotted without and with transverse shear. This time, the flexural stiffness matrix of the sandwich plate is considered totally isotropic ($D_{\alpha\beta}=D$), but we still use the transverse shear stiffness $S$ of the sandwich. The Warburton coefficients used are the special formulations developed for the case of a square plate. The coefficients are more accurate for particular combinations of rows. There is always no difference between Warburton and Blevins computations. Two numerical models are also developed: 2D and 3D. In the 3D model, the honeycomb layer is replaced by a CHEXA8 element, called a PSOLID. There is no difference between the results of the two numerical models. Once again, the analytical and numerical results are still in good agreement. For this plate, an experimental modal analysis has been carried out. The excitation is at the corner point 1 (see Fig. 17). The values of the modal frequencies coming from the experimental analysis are localized with triangular points and marked in italics. Analytical and numerical frequencies are close to these few experimental frequencies. The isotropic square plate is a complex application and is even with composite materials. There are some differences between the analytical rows and the observed numerical or experimental rows. For example, the particular analytical mode (0/4+0/4) which is theoretically the addition of the two simple deformations of modes (0-4) and (4-0) appears in simple form (4-0) for the experiments and the numerical analysis.

The general Eq. (83) and the above applications for free conditions (rectangular and square plates) show that if we have frequencies with transverse shear and $C_{\varphi,m,n}^S$, we can by knowing $S$ and identifying the rows $(m,n)$, return to the frequencies without transverse shear:

$$f_{\varphi,m,n}^{\text{without transverse shear } S} = \frac{\pi}{2a^2} \left( \frac{\langle \lambda_{\varphi,m,n}^{2_{\text{Warburton}}}; D_{\alpha\beta} \rangle}{m_s} \right)^{1/2} = \frac{f_{\varphi,m,n}^S}{C_{\varphi,m,n}^S} \quad (85)$$

The square of the modal frequency without transverse shear represents a linear combination with the elements of the bending matrix:

$$f_{\varphi,m,n}^2 = a_0 \left[ a_1^{\varphi,m,n}, a_2^{\varphi,m,n}, a_3^{\varphi,m,n}, a_4^{\varphi,m,n} \right] \cdot \begin{pmatrix} D_{11} \\ D_{22} \\ H \\ D_{66} \end{pmatrix}, \quad (86)$$

with

$$a_0 = \frac{\pi^2}{4m_s}, \quad a_1^{\varphi,m,n} = \frac{G_{x_{\varphi,m}}^4}{a^4}, \quad a_2^{\varphi,m,n} = \frac{G_{y_{\varphi,n}}^4}{b^4}, \quad a_3^{\varphi,m,n} = \frac{2}{a^2 b^2} H_{x_{\varphi,m}} H_{y_{\varphi,n}}, \quad a_3^{\varphi,m,n} = \frac{4}{a^2 b^2} (j_{x_{\varphi,m}} j_{y_{\varphi,n}} - H_{x_{\varphi,m}} H_{y_{\varphi,n}}). \quad (87)$$
<table><caption>**Table 5**

Warburton coefficients for a plate under free boundary conditions.</caption><thead><tr><th>m</th><th>G<sub>xFP,m</sub></th><th>H<sub>xFP,m</sub></th><th>J<sub>xFP,m</sub></th><th>n</th><th>G<sub>yFP,m</sub></th><th>Hy<sub>FP,m</sub></th><th>Jy<sub>FP,m</sub></th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1</td><td>0</td><td>0</td><td>12/π<sup>2</sup></td><td>1</td><td>0</td><td>0</td><td>12/π<sup>2</sup></td></tr><tr><td>2</td><td>1.506</td><td>1.248</td><td>5.017</td><td>2</td><td>1.506</td><td>1.248</td><td>5.017</td></tr><tr><td>3,4,5,...</td><td>m<sup>-1/2</sup></td><td>(m<sup>-1/2</sup>)<sup>2</sup> [1 - <sup>2</sup>&frasl;<sub>(m<sup>-1/2</sup>)π</sub>]</td><td>(m<sup>-1/2</sup>)<sup>2</sup> [1 + <sup>6</sup>&frasl;<sub>(m<sup>-1/2</sup>)π</sub>]</td><td>3,4,5,...</td><td>n<sup>-1/2</sup></td><td>(n<sup>-1/2</sup>)<sup>2</sup> [1 - <sup>2</sup>&frasl;<sub>(n<sup>-1/2</sup>)π</sub>]</td><td>(n<sup>-1/2</sup>)<sup>2</sup> [1 + <sup>6</sup>&frasl;<sub>(n<sup>-1/2</sup>)π</sub>]</td></tr></tbody></table>

To solve $D_{11}$, $D_{22}$, H and $D_{66}$, from Eq. (86), four different frequencies without out-of-plane shear are needed to form a linear system [$f^2$]=cte.[A].[D]. After the inversion of the linear system, we have:

$$ \begin{Bmatrix} D_{11} \\ D_{22} \\ H \\ D_{66} \end{Bmatrix} = \frac{1}{a_0} \cdot \begin{bmatrix} a_1^{\phi,m_1,n_1} & a_2^{\phi,m_1,n_1} & a_3^{\phi,m_1,n_1} & a_4^{\phi,m_1,n_1} \\ a_1^{\phi,m_2,n_2} & a_2^{\phi,m_2,n_2} & a_3^{\phi,m_2,n_2} & a_4^{\phi,m_2,n_2} \\ a_1^{\phi,m_3,n_3} & a_2^{\phi,m_3,n_3} & a_3^{\phi,m_3,n_3} & a_4^{\phi,m_3,n_3} \\ a_1^{\phi,m_4,n_4} & a_2^{\phi,m_4,n_4} & a_3^{\phi,m_4,n_4} & a_4^{\phi,m_4,n_4} \end{bmatrix}^{-1} \cdot \begin{Bmatrix} f_{m_1,n_1}^2 \\ f_{m_2,n_2}^2 \\ f_{m_3,n_3}^2 \\ f_{m_4,n_4}^2 \end{Bmatrix} \quad (88) $$

The frequencies should be chosen such that the above matrix [A] is not singular. The last component $D_{12}$ is determined with $D_{12}=H-2D_{66}$.

## 6.2. Rectangular plate under the free boundary conditions

It is also interesting to look in detail at the Warburton coefficients for a rectangular plate under free boundary conditions. Table 5 gives the coefficients as a function of the rows m and n.

Under these free conditions, the modal frequencies with transverse shear are expressed as follows:

$$ (f_{FP,m,n}^S)^2 = \frac{a_0 (a_1^{FP,m,n} D_{11} + a_2^{FP,m,n} D_{22} + a_3^{FP,m,n} H + a_4^{FP,m,n} D_{66}))}{1 + (\pi^2/S) (D_{11} (m/a)^2 + D_{22} (n/b)^2)} \quad (89) $$

$D_{11}$, $D_{22}$, $D_{12}$ and $D_{66}$ can be determined from the previous method using four different frequencies, but only after we have corrected the effect of transverse shear. The correction $C_{FP,m,n}^S$ requires the knowledge of S but also $D_{11}$ and $D_{22}$. The elements $D_{11}$ and $D_{22}$ can be determined more directly with the modes called beam modes (m,0) and (0,n), from the frequencies with transverse shear.

### 6.2.1. Beam modes in the x direction

If we focus on the so-called beam modes in the x direction, or the modes of type (m,0), in this case $G_{yFP,0} = H_{yFP,0} = J_{yFP,0} = 0$. Thus, we observe that $a_2^{FP,m,0} = a_3^{FP,m,0} = a_4^{FP,m,0} = 0$ and that it remains just one term $a_1^{FP,m,0} = G_{xFP,m}^4/a^4$. The frequencies of beam modes in the x direction are reduced to:

$$ (f_{FP,m,0}^S)^2 = \frac{a_0 a_1^{FP,m,0} D_{11}}{1 + (\pi^2/S) (D_{11} (m/a)^2)}, \quad (90) $$

from which $D_{11}$ can be deduced:

$$ D_{11} = \frac{(f_{FP,m,0}^S)^2}{a_0 a_1^{FP,m,0} - (\pi^2/S)(m/a)^2 (f_{FP,m,0}^S)^2}, \quad (91) $$

Note: If $m=1$, $a_1^{FP,1,0} = 0$ and $f_{FP,1,0}^S = 0$, it is a rigid body mode and in this case it does not work. It is therefore imperative to take $n=0$ and $m > 1$.

### 6.2.2. Beam modes in the y direction

In the same manner, we find $G_{xFP,0} = H_{xFP,0} = J_{xFP,0} = 0$, $a_1^{FP,0,n} = a_3^{FP,0,n} = a_4^{FP,0,n} = 0$ and with $a_2^{FP,0,n} = G_{yFP,n}^4/b^4$:

$$ (f_{FP,0,n}^S)^2 = \frac{a_0 a_2^{FP,0,n} D_{22}}{1 + (\pi^2/S) (D_{22} (n/b)^2)}, \quad (92) $$
and

$$D_{22} = \frac{\left(f_{FP,0,n}^{S}\right)^{2}}{a_{0} a_{2}^{FP,0,n} - \left(\pi^{2}/S\right)\left(n/b\right)^{2}\left(f_{FP,0,n}^{S}\right)^{2}} \quad (93)$$

Formulas, from Eqs. (91) and (93), then from Eq. (88), have been tested for analytical frequencies of Plate A. They give the original values of $D_{11}$, $D_{22}$, $D_{12}$ and $D_{66}$. From the determination of all elements of the bending matrix, the physical properties of the unidirectional composite ($E_{1}^{UD}$, $E_{2}^{UD}$, $G_{12}^{UD}$ and $\nu_{12}^{UD}$) are now available (see at the end of Section 4.2 and the application in Appendix A).

## 7. Conclusions

The aim of the present work is to propose an analytical formula of eigenmodes for rectangular honeycomb sandwich panels that takes into account the transverse shear of the plate under two boundary conditions all around the plate (simply supported or free). The formulation that has been developed is however only valid for thin or moderately thick plates due to the use of Mindlin's theory. The order of the differential equation being also reduced from 6th to 4th, the simplified model is only valid for panels in which the face sheet thickness is very thin in comparison to the core thickness. These developments are mainly for the low and medium frequency ranges. In addition, the out-of-plane shear stiffness of the honeycomb core should be low, which as well causes a low out-of-plane shear stiffness S of the panel. The bending stiffness of entire sandwich $D_{\alpha\beta}$ should be high and mass per unit area low. The model proposed can fail to predict a finite critical frequency for a sandwich panel type. However, as discussed in Section 4.5, an improved model based on a 6th order differential equation governing the vibration of the structure would ensure a finite critical frequency.

In the structural dynamic domain, the behavior of a symmetric honeycomb sandwich with faces made of laminate is sometimes similar to that of a honeycomb sandwich with isotropic skins. In both cases, the out-of-plane shear stiffness is mainly due to the different thicknesses of the sandwich (faces, honeycomb) and the only transverse shear moduli of the honeycomb. If the laminated skins are thin compared to the honeycomb, which is usually the case, the global transverse shear stiffness can be estimated from the geometric transverse shear stiffness of the honeycomb. In the frequency domain, the global transverse shear stiffness acts on other physical quantities: the module of the wavenumber, the modal frequencies, the modal density and the mechanical impedance. A transition frequency appears at a low frequency. It separates the domain of the classical plate without transverse shear from that with transverse shear. Numerical results presented in this paper are in agreement with the analytical and experimental results.

We also show that the mass per unit area can easily change the value of the theoretical critical frequency. This is particularly important and delicate when panels become lighter. Locally concentrated additional masses (such as glue or densification of the contour) increase the weight. They disrupt the true estimate of the mass per unit area, acting on the acoustic behavior of the panel. That is why it is sometimes difficult to compare acoustic experiments with theoretical results. The choice of the model used to represent the differential equation displacement is also important.

The modal frequencies with transverse shear stiffness for a rectangular sandwich panel under simply supported boundary conditions are well known, but under free boundary conditions it is more difficult. For experiments however, these conditions are the most widely used methods. After analytical and numerical works, we have developed a simple formula suggested in this article for the modal frequencies. Finally, with all parameters controlled, it is possible to obtain some honeycomb and composite face properties from dynamic measurement and through a special process. From the experimental point of view, the transverse shear stiffness S can be estimated by simply measuring the mechanical impedance or the modal density, or the modulus of the wavenumber, which is more complicated. We can trace back to the geometric transverse shear modulus $G_{g}^{H}$ of the honeycomb. With the identification of frequencies from a modal analysis of the sandwich plate and then the determination of all elements of the bending matrix $D_{\alpha\beta}$, the physical properties of the unidirectional composite ($E_{1}^{UD}$, $E_{2}^{UD}$, $G_{12}^{UD}$ and $\nu_{12}^{UD}$) are available. This is possible provided that Young's moduli of laminated skins are large enough compared to those of the honeycomb. It is like having low bending stiffness of the core in front of those produced by the two separated face sheets in the sandwich.

## Appendix A. Application to identify the mechanical properties of the UD from the bending stiffness of the total honeycomb sandwich

With our own laminate [UD(0°,45°, 90°, -45°)$_{S}$, Honeycomb]$_{S}$

### A.1. Bending stiffness with the full equation

The bending stiffness of the honeycomb sandwich can be computed with the honeycomb (Subscript H) and unidirectional composite (Subscript UD) properties, using Eq. (A.1):

$$\{D_{\alpha\beta}\}_{SHC} = \frac{2}{3} [\mathbf{C}_1 \cdot \{Q_{\alpha\beta}^*\}_{UD} + d_0 \cdot \mathbf{C}_0 \cdot \{Q_{\alpha\beta}^*\}_{H}] \quad (A.1)$$
<table><caption>Table A1</caption><thead><tr><th>{D<sub>ij</sub>} = (2/3)</th><th></th><th>C<sub>1</sub></th><th></th><th>{Q<sup>*</sup><sub>ij</sub>}UD</th><th>+</th><th>(2/3)·d<sub>ij</sub></th><th></th><th>C<sub>0</sub></th><th></th><th>{Q<sup>*</sup><sub>ij</sub>}U</th></tr></thead><tbody><tr><td>D<sub>11</sub></td><td></td><td>4.697E-08</td><td>4.685E-08</td><td>3.124E-08</td><td>6.248E-08</td><td>1.343E+11</td><td>1</td><td>0</td><td>0</td><td>1.421E+06</td><td>4621.29</td></tr><tr><td>D<sub>22</sub></td><td></td><td>4.685E-08</td><td>4.697E-08</td><td>3.124E-08</td><td>6.248E-08</td><td>7.738E+09</td><td>0</td><td>1</td><td>0</td><td>1.421E+06</td><td>4611.40</td></tr><tr><td>D<sub>12</sub> = 6.667E-01</td><td>×</td><td>1.562E-08</td><td>1.562E-08</td><td>9.382E-08</td><td>-6.248E-08</td><td>2.244E+09</td><td>+</td><td>1.404E-07</td><td>×</td><td>1.279E+06</td><td>= 1490.01</td></tr><tr><td>D<sub>66</sub></td><td></td><td>1.562E-08</td><td>1.562E-08</td><td>-3.124E-08</td><td>6.258E-08</td><td>3.100E+09</td><td></td><td></td><td></td><td>1.620E+05</td><td>1561.23</td></tr><tr><td>D<sub>16</sub></td><td></td><td>-1.76E-11</td><td>1.76E-11</td><td>-1.91E-24</td><td>-3.81E-24</td><td></td><td>0</td><td>0</td><td>0</td><td>1</td><td>-1.48</td><td>-1.48</td></tr><tr><td>D<sub>26</sub></td><td></td><td>-1.76E-11</td><td>1.76E-11</td><td>1.91E-24</td><td>-3.81E-24</td><td></td><td>0</td><td>0</td><td>0</td><td>0</td><td>-1.48</td><td>-1.48</td></tr></tbody></table>

<table><caption>Table A2</caption><thead><tr><th>{D<sub>ij</sub>}reduced</th><th>=</th><th>(2/3)</th><th></th><th>C<sub>1</sub><sup>reduced</sup></th><th></th><th></th><th>{Q<sup>*</sup><sub>ij</sub>}UD</th><th></th><th></th><th>{Q<sup>*</sup><sub>ij</sub>}U</th></tr></thead><tbody><tr><td>D<sub>11</sub></td><td></td><td></td><td></td><td>4.697E-08</td><td>4.685E-08</td><td>3.124E-08</td><td>6.248E-08</td><td>1.343E+11</td><td></td><td>4621.09</td></tr><tr><td>D<sub>22</sub></td><td>=</td><td>6.667E-01</td><td>×</td><td>4.685E-08</td><td>4.697E-08</td><td>3.124E-08</td><td>6.248E-08</td><td>7.738E+09</td><td>=</td><td>4611.20</td></tr><tr><td>D<sub>12</sub></td><td></td><td></td><td></td><td>1.562E-08</td><td>1.562E-08</td><td>9.382E-08</td><td>-6.248E-08</td><td>2.244E+09</td><td></td><td>1489.83</td></tr><tr><td>D<sub>66</sub></td><td></td><td></td><td></td><td>1.562E-08</td><td>1.562E-08</td><td>-3.124E-08</td><td>6.258E-08</td><td>3.100E+09</td><td></td><td>1561.20</td></tr></tbody></table>

<table><caption>Table A3</caption><thead><tr><th>{Q<sup>*</sup><sub>ij</sub>}UD</th><th>=</th><th>(3/2)</th><th></th><th>C<sub>1</sub><sup>reduced-1</sup></th><th></th><th></th><th>{D<sub>ij</sub>}reduced</th><th></th><th></th><th>{Q<sup>*</sup><sub>ij</sub>}U</th></tr></thead><tbody><tr><td>Q<sup>*</sup><sub>11</sub> UD</td><td></td><td></td><td></td><td>5.603E+09</td><td>-2.930E+09</td><td>-2.665E+09</td><td>-5.329E+09</td><td>4621.09</td><td></td><td>1.343E+11</td></tr><tr><td>Q<sup>*</sup><sub>22</sub> UD</td><td>=</td><td>1.5</td><td>×</td><td>-2.930E+09</td><td>5.603E+09</td><td>-2.665E+09</td><td>-5.329E+09</td><td>4611.20</td><td>=</td><td>7.738E+09</td></tr><tr><td>Q<sup>*</sup><sub>12</sub> UD</td><td></td><td></td><td></td><td>-1.332E+09</td><td>-1.332E+09</td><td>2.673E+09</td><td>5.329E+09</td><td>1489.83</td><td></td><td>2.244E+09</td></tr><tr><td>Q<sup>*</sup><sub>66</sub> UD</td><td></td><td></td><td></td><td>-1.332E+09</td><td>-1.332E+09</td><td>2.665E+09</td><td>5.337E+09</td><td>1561.20</td><td></td><td>3.100E+09</td></tr></tbody></table>
<table><caption>**Table A4**</caption><thead><tr><th></th><th>Q<sup>UD</sup></th><th></th><th></th><th>S<sup>UD</sup>=Q<sup>UD</sup><sup>-1</sup></th><th></th></tr></thead><tbody><tr><td>1.343E+11</td><td>2.244E+09</td><td>0.000E+00</td><td>7.485E-12</td><td>-2.171E-12</td><td>0.000E+00</td></tr><tr><td>2.244E+09</td><td>7.738E+09</td><td>0.000E+00</td><td>-2.171E-12</td><td>1.299E-10</td><td>0.000E+00</td></tr><tr><td>0.000E+00</td><td>0.000E+00</td><td>3.100E+09</td><td>0.000E+00</td><td>0.000E+00</td><td>3.226E-10</td></tr></tbody></table>

<table><caption>**Table A5**</caption><tbody><tr><td>E<sup>UD</sup><sub>1UD</sub> = 1/S<sup>UD</sup><sub>1B</sub> =</td><td>1.336E+11</td><td>Pa</td><td>Whether</td><td>133.6</td><td>GPa</td></tr><tr><td>E<sup>UD</sup><sub>2UD</sub> = 1/S<sup>UD</sup><sub>2B</sub> =</td><td>7.700E+09</td><td>Pa</td><td>Whether</td><td>7.7</td><td>GPa</td></tr><tr><td>G<sup>UD</sup><sub>1B</sub> = 1/S<sup>UD</sup><sub>1B</sub> =<br/>v<sup>UD</sup><sub>1B</sub> = -S<sup>UD</sup><sub>1B</sub>/S<sup>UD</sup><sub>1B</sub> =</td><td>3.100E+09</td><td>Pa</td><td>Whether</td><td>3.1</td><td>GPa</td></tr><tr><td>v<sup>UD</sup><sub>21</sub> = -S<sup>UD</sup><sub>21</sub>/S<sup>UD</sup><sub>22</sub> =</td><td>0.29</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>0.02</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Matrix **C**<sub>1</sub> is computed with the knowledge of the parameters *d*<sub>r</sub> and *θ*<sub>r</sub> associated with each ply *r*. They depend only on the manufacture specifications (honeycomb and UD ply thicknesses, UD orientations). This is a purely geometrical computation (Table A1).

With the full Eq. (A.1), we have $D_{11}=4621$ N m, $D_{22}=4611$ N m, $D_{12}=1490$ N m, $D_{66}=1561$ N m and $D_{16}, D_{26} \approx 0$.

## A.2. Bending stiffness with the simplified equation

If the honeycomb material properties do not matter, Eq. (A.1) can be simplified to Eq. (A.2a).

$$ \{D_{\alpha\beta}\}_{SHC}^{\text{simplified}} \approx \frac{2}{3}\mathbf{C}_1 \cdot \{Q_{\alpha\beta}^*\}_{UD} \qquad (\text{A.2a}) $$

Since $D_{16}$ and $D_{26}$ are close to zero, {$D_{\alpha\beta}$} and **C**<sub>1</sub> matrices are voluntary reduced to the four first lines (see Table A2).

$$ \{D_{\alpha\beta}\}_{SHC}^{\text{reduced}} \approx \frac{2}{3}\mathbf{C}_1^{\text{reduced}} \cdot \{Q_{\alpha\beta}^*\}_{UD} \qquad (\text{A.2b}) $$

With the reduced equation, we find $D_{11}=4621$ N m, $D_{22}=4611$ N m, $D_{12}=1490$ N m and $D_{66}=1561$ N m. There is a slight difference compared to the previous results.

## A.3. Identification of the mechanical properties of the UD

If $D_{11}, D_{22}, D_{12}$ and $D_{66}$ are known, for example from an experimental test, it is possible to identify the properties of the UD. The process requires {$Q_{\alpha\beta}^*$}_{UD} to be solved from the reduced Eq. (A.2b) and then the calculation of the flexibility matrix $S_{UD}^*$ equal to the inverse of $Q_{UD}^*$ matrix. At the end of the process, the properties of the UD are extracted from the $S_{UD}^*$ elements (Tables A3–A5).

First step (Table A3)

$$ \{Q_{\alpha\beta}^*\}_{UD} = \frac{3}{2}\mathbf{C}_1^{\text{reduced} - 1} \cdot \{D_{\alpha\beta}\}_{SHC}^{\text{reduced}} \qquad (\text{A.3a}) $$

Second step (Table A4)

$$ S_{UD}^* = Q_{UD}^* - 1 \qquad (\text{A.3b}) $$

Third and last step (Table A5)

With the knowledge of $D_{11}, D_{22}, D_{12}, D_{66}$ and matrix **C**<sub>1</sub>, the mechanical properties of the UD are identified.

## References

[1] J.P.D. Wilkinson, Modal densities of certain shallow structural elements, *Journal of the Acoustical Society of America* 43 (2) (1968) 245-251.

[2] R.H. Lyon, R.G. DeJong, *Theory and Application of Statistical Energy Analysis*, second ed. Butterworth-Heinemann, Newton, 1995.

[3] L.L. Erickson, Modal densities of sandwich panels: theory and experiment, *Shock & Vibration Bulletin* 39 (3) (1969) 1-16.

[4] D.J. Mead, The damping properties of elastically supported sandwich plates, *Journal of Sound and Vibration* 24 (3) (1972) 275-295.

[5] B.L. Clarkson, M.F. Ranky, Modal density of honeycomb plates, *Journal of Sound and Vibration* 91 (1) (1983) 103-118.

[6] K. Renji, P.S. Nair, S. Narayanan, Modal density of composite honeycomb sandwich panels, *Journal of Sound and Vibration* 195 (5) (1996) 687-699.

[7] R.D. Blevins, *Formulas for Natural Frequency and Mode Shape*, Editor Van Nosirand Reinhold or Krieger Company, New York, 1979.

[8] T. Lauwagie, H. Sol, W. Heylen, G. Roebben, Determination of in-plane elastic properties of the different layers of laminated plates by means of vibration testing and model updating, *Journal of Sound and Vibration* 274 (2004) 529-546.

[9] S.D. Yu, W.L. Cleghorn, Free flexural vibration analysis of symmetric honeycomb panels, *Journal of Sound and Vibration* 284 (2005) 189-204.

[10] G. Kurtze, B. Watters, New wall design for high transmission loss or high damping, *Journal of the Acoustical Society of America* 31 (6) (1959) 739-748.

[11] D. Gay, *Matériaux Composites (Chapitre 18: Applications Poutre Sandwich)*, 3rd ed. Edition Hermès, 1997.

[12] J.-M. Berthelot, *Les matériaux Composites – Comportement Mécanique et Analyse des Structures*, 5th ed. Tec&Doc Edition, Paris, 2012.
[13] J.-M. Berthelot, *Composite Materials, Mechanical behavior and Structural Analysis*, Springer Edition, New York, 1999.

[14] G.B. Warburton, The vibration of rectangular plates, *Proc. Inst. Mech. Eng. 168* (12) (1954) 371–381.

[15] C. Lesueur, *Rayonnement acoustique des structures*, Eyrolles Edition, Paris, 1988.

[16] N.G. Stephen, Mindlin plate theory: best shear coefficient and higher spectra validity, *Journal of Sound and Vibration* 202 (4) (1997) 539–553.

[17] N.G. Stephen, On a check on the accuracy of Timoshenko's beam theory, *Journal of Sound and Vibration* 257 (4) (2002) 809–812.

[18] J.R. Hutchinson, Shear coefficients for Timoshenko beam theory, *Journal of Applied Mechanics* 68 (2001) 87–92.

[19] K. Renji, P.S. Nair, S. Narayanan, Critical and coincidence frequencies of flat panels, *Journal of Sound and Vibration* 205 (1) (1997) 19–32.

[20] L.J. Gibson, M.F. Ashby, *Cellular Solids – Structure and Properties* (re-issued 1997), Cambridge University Press, 1988.

[21] X. Carniel, L. Gavric, Study formulations vibration intensity for bending plates, structural intensity and vibrational energy flow, in: *3rd International Congress on Intensity Techniques*, Cetim Senlis, 1990.