**1. Introduction**

Nowadays, along with a strong development of science and technology, there are many new advanced materials appeared, for instance, composite materials, functionally graded materials (FGM), piezoelectric materials, and so on. The studying on dynamic responses of these new materials has been reached grea<sup>t</sup> achievements and attracted numerous scientists all over the world. Moreover, the idea of merging these different materials is considered by engineers to make new structures in order to have specific purposes. For example, the combining of a concrete structure and a steel structure has a lighter weight than a normal concrete structure. Hence, these new types of structures are applied extensively in civil techniques, aerospace, and army vehicles. In this structure, the connecting stub is attached to contact different layers in order to create the compatibility of the horizontal displacement among layers, and it plays an important role in working process of the structure.

For multilayered beams, recently, the Newark's model [1] is considered by many experts such as He et al. [2], Xu and Wang [3,4]. They took into account the shear strain when calculating by using Timoshenko beam theory. Nguyen [5] studied the linear dynamic problems. Silva et al. [6], Schnabl et al. [7] and Nguyen and co-workers [8,9] employed the finite element method (FEM) and analytical method in order to examine linear static analysis of multilayered beams. Huang [10], and Shen [11] studied the linear dynamic response, too. For the nonlinear free vibration can be seen in [12] of Arvin and Bakhtiari-Nejad.

In addition to the Timoshenko beam theory (TBT), the higher-order beam theory (HBT) is also considered, in which The dynamic problem is carried out by Chakrabarti in [13] with FEM. Chakrabarti and colleagues [14,15] analyzed a static problem for two-layer composite beams. The higher-order beam theory (HBT) overcomes a part of the e ffect due to the shear locking coe fficient caused. Otherwise, Subramanian [16] constructed an element based on a displacement field to study the free vibration of the multilayered beam. Li et al. [17] conducted a free vibration analysis by employing the hyperbolic shear deformation theory. Vo and Thai [18] studied static multilayered beams with the improved higher-order beam theory of Shimpi.

In general, most higher-order beam theories (HBT), including higher-order beam theory of Reddy tend to neglect the horizontal deformation of multilayered beams. According to the Kant's opinion, the horizontal stress of sub-layer is caused by the pressure can reduce the dimension from multilayered beam model to the plane stress model. To obtain this thing, Kant [19,20] employed both the higher-order beam theory (HBT) and the horizontal displacement theory by considering approximate displacements in two ways. Thus, he established the mixed two-layer beam with sub-layers, which abides by the higher-order beam theory of Kant proposed by the weak form for the buckling analysis.

A three-dimensional fracture plasticity based on finite element model (FEM) are developed by Yan and coworkers [21] to carry out the ultimate strength respones of SCS sandwich structure under concentrated loads. The static behaviors of beams with di fferent types of cross-section, such as square, C-shaped, and bridge-like sections, were investigated in Carrena's study [22] by assuming that the displacement field is expanded in terms of generic functions, which is the Unified Formulation by Carrera (CUF) [23]. Similarly to mentioned methods, Cinefra et al. [24] used MITC9 shell elements to explore the mechanical behavior of laminated composite plates and shells. Muresan and coworkers [25] examined the study on the stability of thin walled prismatic bars based on the Generalized Beam Theory (GBT), which is an e fficient approach developed by Schardt [26]. Yu et al. [27] employed the Variational Asymptotic Beam Section Analysis (VABS) for mechanical behavior of various cross-sections such as elliptic and triangular sections. In [28], we used first-order shear deformation theory to analysis of triple-layer composite plates with layers connected by shear connectors subjected to moving load. Ansari Sadrabadi et al. [29] used analytical methods to investigate a thick-walled cylindrical tube made of a functionally graded material (FGM) and undergoing thermomechanical loads.

For multilayered plate and shell composite structures, there have been many published papers, including static problems, dynamic problems, linear, and nonlinear problems, and so on. However, for the multilayered structure with shear connectors, there are not many papers yet. Based on above mentioned papers, the authors are about to construct the relations of mechanical behavior and the oscillator equation of the multilayered shell. We also study several geometrical and physical parameters, the loading, etc., which e ffect on the vibration of the shell.

The body of this paper is divided into five main sections. Section 1 is the general introduction. We present finite formulations of free vibration and forced vibration analysis of three-layer composite shell with shear connectors in Section 2. The numerical results of vibration and forced vibrations are discussed in Sections 3 and 4. Section 5 gives some major conclusions.

### **2. Finite Element Formulations**

### *2.1. Equation of Motion of the Shell Element*

Consider a three-layer composite shell with shear connectors as shown in Figure 1.

**Figure 1.** The model of three-layer composite shell with shear connectors, (**a**) shell model with shear connectors, and (**b**) finite element model.

The composite shell consists of three layers, including the top layer (*t*), the bottom layer (*b*) and the middle layer (*c*); these layers are connected with one another by shear connectors, and they can be made of the same materials or different materials. These three layers can slide relatively with one another at the contact surfaces, and there is no delamination phenomenon at all. All three layers of the shell are set in the local coordinates *Oxyt*, *Oxyc*, and *Oxyb*, respectively. The total thickness of the shell is divided into six small part *h*1, *h*2, *h*3, *h*4, *h*5, *h*6 as shown in Figure 1; *ut*0, *uc*0, and *ub*0 represent displacements in *x* direction; *vb*0 represents the displacement in y direction at the neutral surface of each layer.

According to Mindlin plate theory, displacements *u, v, w* at a point (*xk, yk, zk*) of layer *k* are expressed as follows:

$$\begin{cases} u\_k = u\_{k0}(\mathbf{x}\_{k'} y\_k) + z\_k \boldsymbol{\uprho}\_k(\mathbf{x}\_{k'} y\_k) \\ v\_k = v\_{k0}(\mathbf{x}\_{k'} y\_k) + z\_k \boldsymbol{\uprho}\_k(\mathbf{x}\_{k'} y\_k) \\ w\_k = w(\mathbf{x}\_k, y\_k) \end{cases} (k = t, c, b) \tag{1}$$

where ϕ*k* and ψ*k* are the transverse normal rotations of the *xk* and *yk* directions.

The relative movements among the contact surfaces are defined by the following equations For the layer *t* and layer *c* we have

$$\begin{cases} \boldsymbol{u}\_{\rm lc} = \boldsymbol{u}\_{\rm l}(\mathbf{x}\_{\rm l}, \mathbf{y}\_{\rm l}, \mathbf{h}\_2) - \boldsymbol{u}\_{\rm c}(\mathbf{x}\_{\rm c}, \mathbf{y}\_{\rm c} - \mathbf{h}\_3) \\ \boldsymbol{v}\_{\rm lc} = \boldsymbol{v}\_{\rm l}(\mathbf{x}\_{\rm l}, \mathbf{y}\_{\rm l}, \mathbf{h}\_2) - \boldsymbol{v}\_{\rm c}(\mathbf{x}\_{\rm c}, \mathbf{y}\_{\rm c} - \mathbf{h}\_3) \end{cases} \tag{2}$$

And for layer *c* and layer *b* we have:

$$\begin{cases} \boldsymbol{u}\_{cb} = \boldsymbol{u}\_{c}(\boldsymbol{x}\_{c}, \boldsymbol{y}\_{b^{\prime}}, \boldsymbol{h}\_{4}) - \boldsymbol{u}\_{b}(\boldsymbol{x}\_{b^{\prime}}, \boldsymbol{y}\_{b^{\prime}} - \boldsymbol{h}\_{5}) \\ \boldsymbol{v}\_{cb} = \boldsymbol{v}\_{c}(\boldsymbol{x}\_{c^{\prime}}, \boldsymbol{y}\_{c^{\prime}}, \boldsymbol{h}\_{4}) - \boldsymbol{v}\_{b}(\boldsymbol{x}\_{b^{\prime}}, \boldsymbol{y}\_{b^{\prime}} - \boldsymbol{h}\_{5}) \end{cases} \tag{3}$$

Note that at the contact surfaces, we have:

$$\begin{cases} z\_t = h\_2; z\_c = -h\_3\\ z\_t = h\_4; z\_b = -h\_5 \end{cases} \tag{4}$$

with *h*4 = *h*3 = *hc*2 .

> From Equations (1)–(4), we get:

$$\begin{cases} \boldsymbol{u}\_{t\boldsymbol{\varepsilon}} = \boldsymbol{u}\_{t0} - \boldsymbol{u}\_{t0} + h\_2 \boldsymbol{q}\_t + h\_3 \boldsymbol{q}\_{\boldsymbol{\varepsilon}}\\ \boldsymbol{v}\_{t\boldsymbol{\varepsilon}} = \boldsymbol{v}\_{t0} - \boldsymbol{v}\_{t0} + h\_2 \boldsymbol{\psi}\_t + h\_3 \boldsymbol{\psi}\_{\boldsymbol{\varepsilon}} \end{cases} \tag{5}$$

$$\begin{cases} u\_{cb} = u\_{c0} - u\_{l0} + h\_4 \varphi\_c + h\_5 \varphi\_b \\\ v\_{cb} = v\_{c0} - v\_{l0} + h\_4 \psi\_c + h \mathfrak{g} \psi\_b \end{cases} \tag{6}$$

The relation between strain and displacement of each layer is expressed as follows For the layer *k*, we have:

$$\begin{array}{l} \varepsilon\_{\mathbf{k}\mathbf{r}} = \frac{\partial u\_{\mathbf{k}}}{\partial \mathbf{x}} = \frac{\partial u\_{\mathbf{k}0}}{\partial \mathbf{x}} + \frac{u\_{\mathbf{0}}}{\mathbf{R}\_{\mathbf{r}}} + z\_{\mathbf{k}} \frac{\partial \mathbf{v}\_{\mathbf{k}}}{\partial \mathbf{x}};\\ \varepsilon\_{\mathbf{k}\mathbf{y}} = \frac{\partial v\_{\mathbf{0}}}{\partial \mathbf{y}} = \frac{\partial v\_{\mathbf{0}0}}{\partial \mathbf{y}} + \frac{u\_{\mathbf{0}}}{\mathbf{R}\_{\mathbf{y}}} + z\_{\mathbf{k}} \frac{\partial \psi\_{\mathbf{k}}}{\partial \mathbf{y}};\\ \gamma\_{\mathbf{k}\mathbf{xy}} = \frac{\partial v\_{\mathbf{k}}}{\partial \mathbf{x}} + \frac{\partial u\_{\mathbf{k}}}{\partial \mathbf{y}} = \frac{\partial u\_{\mathbf{0}0}}{\partial \mathbf{y}} + \frac{\partial v\_{\mathbf{0}0}}{\partial \mathbf{x}} + \frac{2u\_{\mathbf{0}}}{\mathbf{R}\_{\mathbf{x}}} + z\_{\mathbf{k}} \Big(\frac{\partial p\_{\mathbf{k}}}{\partial \mathbf{y}} + \frac{\partial \psi\_{\mathbf{k}}}{\partial \mathbf{x}} + \frac{1}{2} \Big(\frac{1}{\mathbf{R}\_{\mathbf{y}}} - \frac{1}{\mathbf{R}\_{\mathbf{x}}} \Big) \Big(\frac{\partial v\_{\mathbf{0}0}}{\partial \mathbf{x}} - \frac{\partial u\_{\mathbf{0}}}{\partial \mathbf{y}} \Big);\\ \gamma\_{\mathbf{k}\mathbf{xz}} = \frac{\partial u\_{\mathbf{0}}}{\partial \mathbf{x}} + \frac{\partial u\_{\mathbf{k}}}{\partial \mathbf{z}\_{\mathbf{k}}} = \frac{\partial u\_{\mathbf{0}}}{\partial \mathbf$$

We can rewrite in a matrix form as follow

$$\mathfrak{e}\_{k} = \left\{ \begin{array}{c} \varepsilon\_{kx} \\ \varepsilon\_{ky} \\ \gamma\_{kxy} \end{array} \right\} = \mathfrak{e}\_{k}^{0} + z\_{k} \mathfrak{e}\_{k}; \ \mathfrak{y}\_{k} = \left\{ \begin{array}{c} \mathcal{V}\_{kyz} \\ \mathcal{V}\_{kxx} \end{array} \right\} \tag{8}$$

in which

$$\begin{aligned} \boldsymbol{\varepsilon}\_{k}^{0} &= \begin{pmatrix} \varepsilon\_{0}^{0} \\ \varepsilon\_{ky}^{0} \\ \gamma\_{ky}^{0} \end{pmatrix} = \begin{pmatrix} \frac{\partial \mathbf{u}\_{0}}{\partial \mathbf{x}} + \frac{\mathbf{u}\_{0}}{\mathbf{R}\_{k}} \\ \frac{\partial \mathbf{v}\_{0}}{\partial \mathbf{y}} + \frac{\mathbf{u}\_{0}}{\mathbf{R}\_{y}} \\ \left(\frac{\partial \mathbf{u}\_{0}}{\partial \mathbf{y}} + \frac{\partial \mathbf{v}\_{0}}{\partial \mathbf{x}}\right) + \frac{2\mathbf{u}\_{0}}{\mathbf{R}\_{ky}} \end{pmatrix}; \mathbf{w}\_{k} &= \begin{pmatrix} \kappa\_{kx} \\ \kappa\_{ky} \\ \kappa\_{ky} \end{pmatrix} = \begin{pmatrix} \frac{\partial \mathbf{v}\_{1}}{\partial \mathbf{x}} \\ \frac{\partial \mathbf{v}\_{2}}{\partial \mathbf{y}} \\ \frac{\partial \mathbf{v}\_{3}}{\partial \mathbf{y}} + \frac{\partial \mathbf{u}\_{3}}{\partial \mathbf{x}} + \mathbf{p}\_{k} \\ \frac{\partial \mathbf{v}\_{4}}{\partial \mathbf{y}} - \frac{\mathbf{u}\_{0}}{\mathbf{R}\_{y}} + \frac{\partial \mathbf{u}\_{4}}{\partial \mathbf{y}} + \mathbf{p}\_{k} \end{pmatrix} = \begin{pmatrix} \frac{\partial \mathbf{v}\_{1}}{\partial \mathbf{x}} \\ \frac{\partial \mathbf{v}\_{2}}{\partial \mathbf{y}} \\ \frac{\partial \mathbf{v}\_{3}}{\partial \mathbf{y}} + \frac{\partial \mathbf{u}\_{3}}{\partial \mathbf{x}} + \mathbf{p}\_{k} \end{pmatrix} \end{aligned} \tag{9}$$

The relation between stress and strain of layer *k* is expressed as followIs necessary bild?

$$
\sigma\_k = \mathbf{D}\_k \boldsymbol{\varepsilon}\_k; \ \mathbf{r}\_k = \frac{5}{6} \mathbf{G}\_k \mathbf{y}\_k \tag{10}
$$

in which **D***k*, **G***k* are the bending rigidity and shear rigidity of layer *k*, respectively, and ν*k* is the Poisson ratio of layer *k*.

$$\mathbf{D}\_{k} = \frac{E\_{k}}{1 - v^{2}} \begin{bmatrix} 1 & \nu\_{k} & 0 \\ \nu\_{k} & 1 & 0 \\ 0 & 0 & (1 - \nu\_{k})/2 \end{bmatrix}; \mathbf{G}\_{k} = \frac{E\_{k}}{2(1 + \nu\_{k})} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \tag{11}$$

In this work, the thickness of the shell is thin or medium (*h* = *a*100 ÷ *a*10 , *a* is short edge), we employ the 8-node isoparametric element, each node has 13 degrees of freedom, three layers have the same displacement in the *z*-direction (Figure 2), the degree of freedom of node *i* is **q***ie* and the total degree of freedom of the shell element **q***e* is defined as follow.

$$\mathbf{q}\_e^i = \begin{pmatrix} \ u\_{t0i} & v\_{t0i} & q\_{ti} & \psi\_{ti} & u\_{c0i} & v\_{c0i} & q\_{ci} & \psi\_{ci} & u\_{l0i} & v\_{l0i} & q\_{li} & \psi\_{li} & w \end{pmatrix}^T; i = 1 \div 8. \tag{12}$$

$$\mathbf{q}\_{\boldsymbol{\epsilon}} = \left\{ \begin{array}{cccccccc} \mathbf{q}\_{\boldsymbol{\epsilon}}^{1} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{2} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{3} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{4} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{5} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{6} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{7} & \mathbf{q}\_{\boldsymbol{\epsilon}}^{8} \\ \dots & & & & \end{array} \right\}^{\mathrm{T}} \tag{13}$$

$$\begin{aligned} \mu\_{k0} &= \sum\_{i=1}^{8} N\_i(\boldsymbol{\xi}, \boldsymbol{\eta}) u\_{k0i}; \; \nu\_{k0} = \sum\_{i=1}^{8} N\_i(\boldsymbol{\xi}, \boldsymbol{\eta}) v\_{k0i} \\ \rho\_k &= \sum\_{i=1}^{8} N\_i(\boldsymbol{\xi}, \boldsymbol{\eta}) \rho\_{ki}; \; \psi\_k = \sum\_{i=1}^{8} N\_i(\boldsymbol{\xi}, \boldsymbol{\eta}) \psi\_{ki}; \; w = \sum\_{i=1}^{8} N\_i(\boldsymbol{\zeta}, \boldsymbol{\eta}) w\_i \end{aligned} \tag{14}$$

in which *Ni* (*i* = 1 ÷ 8) can be defined as in [28].

**Figure 2.** Degrees of freedom of the node in the eight-node shell element.

By substituting in the expression for verifying displacement of element we have:

$$\begin{cases} \mathfrak{e}\_{k} = \left( \mathbf{B}\_{k}^{0} + \boldsymbol{z}\_{k} \mathbf{B}\_{k}^{1} \right) \mathbf{q}\_{x} \\ \mathbf{y}\_{k} = \mathbf{S}\_{k} \mathbf{q}\_{c} \end{cases} \quad (k = \boldsymbol{t}, \boldsymbol{c}, \boldsymbol{b}) \tag{15}$$

in which **B**0*k*; **B**1*k*; **S***k* are defined as follows

$$
\begin{array}{l}
\mathbf{B}\_{k}^{0} = \left[\begin{array}{cccccc}\mathbf{B}\_{k1}^{0} & \mathbf{B}\_{k2}^{0} & \mathbf{B}\_{k3}^{0} & \mathbf{B}\_{k4}^{0} & \mathbf{B}\_{k5}^{0} & \mathbf{B}\_{k6}^{0} & \mathbf{B}\_{k7}^{0} & \mathbf{B}\_{k8}^{0} \\ \mathbf{B}\_{k}^{1} = \left[\begin{array}{cccccc}\mathbf{B}\_{k1}^{1} & \mathbf{B}\_{k2}^{1} & \mathbf{B}\_{k3}^{1} & \mathbf{B}\_{k4}^{1} & \mathbf{B}\_{k5}^{1} & \mathbf{B}\_{k7}^{1} & \mathbf{B}\_{k8}^{1} \\ \mathbf{B}\_{k1}^{1} & \mathbf{B}\_{k2}^{1} & \mathbf{B}\_{k3}^{1} & \mathbf{B}\_{k4}^{1} & \mathbf{B}\_{k5}^{1} & \mathbf{B}\_{k7}^{1} & \mathbf{B}\_{k8}^{1} & \mathbf{B}\_{k9}^{1} \\ \end{array} \right]; \tag{16}
$$

where **<sup>B</sup>**0*ki*, **B**1*ki*and **<sup>S</sup>***ki*can be found in Appendix A

The elastic force of connector stub per unit length is defined by the following equations. For layer *t* and *c* we have:

$$\mathbf{F}\_c^{tc} = \left\{ \begin{array}{c} F\_{cu} \\ F\_{cv} \end{array} \right\}\_{ct} = k\_{tc} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{Bmatrix} u\_{lc} \\ v\_{tc} \end{Bmatrix} = \mathbf{K}\_c^{tc} \mathbf{q}\_c^{tc} \tag{17}$$

With

$$\mathbf{q}\_{\varepsilon}^{\text{lc}} = \left\{ \begin{array}{c} \boldsymbol{\mu}\_{\text{lt}} \\ \boldsymbol{v}\_{\text{lt}} \end{array} \right\} = \left[ \begin{array}{c} \boldsymbol{\mu}\_{\text{l}0} + \boldsymbol{h}\_{2}\boldsymbol{\varrho}\_{\text{l}} - \boldsymbol{u}\_{\text{c}0} + \boldsymbol{h}\_{3}\boldsymbol{\varrho}\_{\text{c}} \\ \boldsymbol{v}\_{\text{l}0} + \boldsymbol{h}\_{2}\boldsymbol{\psi}\_{\text{t}} - \boldsymbol{v}\_{\text{c}0} + \boldsymbol{h}\_{3}\boldsymbol{\psi}\_{\text{c}} \end{array} \right] = \mathbf{N}\_{\text{lt}} \mathbf{q}\_{\varepsilon} = \sum\_{i=1}^{8} (\mathbf{N}\_{\text{l}\varepsilon})\_{i} q\_{\varepsilon}^{i} \tag{18}$$

*Symmetry* **2019**, *11*, 527

in which

$$(\mathbf{N}\_{l\mathbf{t}})\_i = \begin{bmatrix} N\_i & 0 & h\_2 N\_i & 0 & -N\_i & 0 & h\_3 N\_i & 0 & 0 \ 0 \ 0 \ 0 \ 0 \\ 0 & N\_i & 0 & h\_2 N\_i & 0 & -N\_i & 0 & h\_3 N\_i & 0 \ 0 \ 0 \ 0 \ 0 \ 0 \end{bmatrix} \tag{19}$$

For layer *c* and *b* we have

$$\mathbf{F}\_{\varepsilon}^{cb} = \left\{ \begin{array}{c} F\_{c\upsilon} \\ F\_{c\upsilon} \end{array} \right\}\_{cb} = k\_{cb} \left[ \begin{array}{c} 1 & 0 \\ 0 & 1 \end{array} \right] \left\{ \begin{array}{c} u\_{cb} \\ v\_{cb} \end{array} \right\} = \mathbf{K}\_{\varepsilon}^{cb} \mathbf{q}\_{\varepsilon}^{cb} \tag{20}$$

with

$$\mathbf{q}\_{\mathbf{f}}^{cb} = \left\{ \begin{array}{c} \boldsymbol{u}\_{cb} \\ \boldsymbol{v}\_{cb} \end{array} \right\} = \left[ \begin{array}{c} \boldsymbol{u}\_{c0} + \boldsymbol{h}\_{\mathbf{4}} \boldsymbol{q}\_{\mathbf{c}} - \boldsymbol{u}\_{b0} + \boldsymbol{h}\_{\mathbf{5}} \boldsymbol{q}\_{\mathbf{b}} \\ \boldsymbol{v}\_{c0} + \boldsymbol{h}\_{\mathbf{4}} \boldsymbol{\uppsi}\_{\mathbf{c}} - \boldsymbol{v}\_{b0} + \boldsymbol{h}\_{\mathbf{5}} \boldsymbol{\uppsi}\_{\mathbf{b}} \end{array} \right] = \mathbf{N}\_{\mathbf{c}b} \mathbf{q}\_{\mathbf{c}} = \sum\_{i=1}^{8} (\boldsymbol{N}\_{\mathbf{c}b})\_{i} \boldsymbol{q}\_{\mathbf{c}}^{i} \tag{21}$$

in which

$$(\mathbf{N}\_{cb})\_i = \begin{bmatrix} 0 \ 0 \ 0 \ 0 \ N\_i & 0 & h\_4 N\_i & 0 & -N\_i & 0 & h\_5 N\_i & 0 & 0\\ 0 \ 0 \ 0 \ 0 \ 0 & N\_i & 0 & h\_4 N\_i & 0 & -N\_i & 0 & h\_5 N\_i & 0 \end{bmatrix} \tag{22}$$

Here, *ktc* and *kcb* are the shear resistance coefficients of the connector stub per unit length. To obtain the dynamic equation we employ the weak form for each element, we get:

$$\begin{split} \sum\_{k=t,\varepsilon,b\in\mathcal{V}\_{k}} \int \delta\dot{\mathbf{q}}\_{k}^{\mathrm{T}} \rho\_{k}\dot{\mathbf{q}}\_{k}dV\_{k} + \sum\_{k=t,\varepsilon} \int \delta\varepsilon\_{k}^{\mathrm{T}} \sigma\_{k}dV\_{k} + \frac{5}{6} \sum\_{k=t,\varepsilon} \int \delta\gamma\_{k}^{\mathrm{T}} \pi\_{k}dV\_{k} + \sum\_{k=t,\varepsilon\,b\in\mathcal{A}\_{k}} \int \delta\Big(\mathbf{q}\_{t}^{k}\Big)^{\mathrm{T}} \mathbf{F}\_{t}^{k}dA\_{k} \\ - \delta\mathbf{q}\_{t}^{\mathrm{T}} \int \mathbf{N}\_{\mathrm{w}}p(t)dA\_{\mathrm{l}} = 0 \end{split} \tag{23}$$

By substituting Equations (1), (15), (17), and (20) into Equation (23), we obtain the dynamic equation of the shell element as follows:

$$\mathbf{M}\_{\varepsilon}\ddot{\mathbf{q}}\_{\varepsilon} + \mathbf{K}\_{\varepsilon}\mathbf{q}\_{\varepsilon} = \mathbf{F}\_{\varepsilon}(t) \tag{24}$$

with

$$\begin{split} \mathbf{K}\_{c(10\&104)} &= \sum\_{k=c,s,a} \int\_{A\_{k}} \mathbf{B}\_{k}^{0} \mathbf{D}\_{k} \mathbf{B}\_{k}^{0} dA\_{k} + \sum\_{k=c,s,a} \int \left( \mathbf{B}\_{k}^{0} \right)^{\mathrm{T}} \mathbf{D}\_{k1} \mathbf{B}\_{k}^{1} dA\_{k} + \\ &+ \sum\_{k=t,c,b} \int \left( \mathbf{B}\_{k}^{1} \right)^{\mathrm{T}} \mathbf{D}\_{k1} \mathbf{B}\_{k}^{0} dA\_{k} + \sum\_{k=t,c,b} \int \left( \mathbf{B}\_{k}^{1} \right)^{\mathrm{T}} \mathbf{D}\_{k2} \mathbf{B}\_{k}^{1} dA\_{k} + \frac{5}{6} \sum\_{k=t,c,b} \int \mathbf{S}\_{k}^{\mathrm{T}} \mathbf{G}\_{k} \mathbf{S}\_{k} dA\_{k} \\ &+ \int \mathbf{N}\_{tc}^{\mathrm{T}} \mathbf{K}\_{tc}^{c} \mathbf{N}\_{tc} dA\_{tc} + \int \mathbf{N}\_{cb}^{\mathrm{T}} \mathbf{K}\_{cb}^{c} \mathbf{N}\_{cb} dA\_{cb} \end{split} \tag{25}$$

in which

$$(\mathbf{D}\_{k0}; \ \mathbf{D}\_{k1}; \ \mathbf{D}\_{k2}) = \int\_{-h\_k/2}^{h\_k/2} \{1; z\_k; z\_k^2 \} \mathbf{D}\_k \, dz\_k; \mathbf{H}\_k = \int\_{-h\_k/2}^{h\_k/2} \mathbf{G}\_k \, dz\_k \ (k = t\_\prime c, b) \tag{26}$$

$$\mathbf{M}\_{\varepsilon(104 \times 104)} = \sum\_{k=t,\varepsilon,b} \int \int \int^{h\_k/2} \mathbf{L}\_k^\mathrm{T} \rho\_k \mathbf{L}\_k dz\_k dA\_k \tag{27}$$

where **L***k* can be seen in Appendix B

$$\mathbf{F}\_{\mathbf{f}}(t)\_{(104\times1)} = \int\_{A\_{\ell}} p(t) \mathbf{N}\_{\mathbf{w}}^{\mathrm{T}} dA\_{\ell} \tag{28}$$

in which

$$\mathbf{N}\_{\rm w} = \begin{bmatrix} \mathbf{N}\_{\rm w1} & \mathbf{N}\_{\rm w2} & \mathbf{N}\_{\rm w3} & \mathbf{N}\_{\rm w4} & \mathbf{N}\_{\rm w5} & \mathbf{N}\_{\rm w6} & \mathbf{N}\_{\rm w7} & \mathbf{N}\_{\rm w8} \end{bmatrix} \tag{29}$$

with

$$\mathbf{N}\_{\mathbf{w}\dot{j}} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 00000 & N\_{\dot{j}} \end{bmatrix} \tag{30}$$

In the case of taking into account the structural damping, we have the force vibration equation of the shell element as follows:

$$\mathbf{M}\_{\varepsilon}\ddot{\mathbf{q}}\_{\varepsilon} + \mathbf{C}\_{\varepsilon}\dot{\mathbf{q}}\_{\varepsilon} + \mathbf{K}\_{\varepsilon}\mathbf{q}\_{\varepsilon} = \mathbf{F}\_{\varepsilon}(t) \tag{31}$$

in which **C***e* = <sup>α</sup>**M***e* + β**K***e* and α, β are Rayleigh drag coefficients defined in [30,31].

### *2.2. The Di*ff*erential Equation of Vibration*

From the differential equation of vibration of the shell element (Equation (31)), we obtain the differential equation of forced vibration of three-layer composite shell as follows:

$$\mathbf{M}\ddot{\mathbf{q}} + \mathbf{C}\mathbf{q} + \mathbf{K}\mathbf{q} = \mathbf{F}(t) \tag{32}$$

in which **M**, **C**, **K**, **F**(*t*) are the global mass matrix, the global structural damping matrix, the global stiffness matrix and the global load matrix, respectively. These matrices and vectors are assembled from the element matrices and vectors, correspondingly. They are linear differential equations, which have the right-hand side depending on time. In order to solve these equations, we use the Newmark-beta method [31]. The program is coded in the MATLAB (MathWorks, Natick, MA, USA) environment with the following algorithm flowchart of Newmark as shown Figure 3.

**Figure 3.** Algorithm flowchart of Newmark solving the dynamic response problem of the shell.

*Symmetry* **2019**, *11*, 527

> For the free vibration analysis, the natural frequencies can be obtained by solving the equation:

$$\mathbf{M}\ddot{\mathbf{q}} + \mathbf{K}\mathbf{q} = 0\tag{33}$$

or in another form:

$$(\mathbf{K} - \omega^2 \mathbf{M})\mathbf{q} = 0\tag{34}$$

where ω is the natural frequency.

> Flowchart of Newmark-beta method [31] Step1:Determine the firstconditions:

$$\mathbf{q}(0) = \mathbf{q}\_0; \dot{\mathbf{q}}(0) = \dot{\mathbf{q}}\_0 \tag{35}$$

From the first conditions, we obtain:

$$
\ddot{\mathbf{q}}\_0 = \mathbf{M}\_0^{-1} \left( \mathbf{F}\_0 - \mathbf{K}\_0 \mathbf{q}\_0 - \mathbf{C}\_0 \dot{\mathbf{q}}\_0 \right) \tag{36}
$$

Step 2: By approximating **..q***<sup>t</sup>*+Δ*t*, .**q***t*+Δ*t* by **q***<sup>t</sup>*+Δ*t*, we have

$$\begin{aligned} \ddot{\mathbf{q}}\_{t+\Delta t} &= a\_0 (\mathbf{q}\_{t+\Delta t} - \mathbf{q}\_t) - a\_2 \dot{\mathbf{q}}\_t - a\_3 \ddot{\mathbf{q}}\_t \\ \dot{\mathbf{q}}\_{t+\Delta t} &= \dot{\mathbf{q}}\_t + a\_6 \ddot{\mathbf{q}}\_t + a\_7 \ddot{\mathbf{q}}\_{t+\Delta t} \end{aligned} \tag{37}$$

where:

$$\begin{aligned} a\_0 &= \frac{2}{\gamma \Delta t^2}; \ a\_1 = \frac{2\alpha}{\gamma \Delta t}; \ a\_2 = \frac{2}{\gamma \Delta t}; \ a\_3 = \frac{1}{\gamma} - 1; \ a\_4 = \frac{2\alpha}{\gamma} - 1;\\ a\_5 &= \left(\frac{\alpha}{\gamma} - 1\right) \Delta t; \ a\_6 = (1 - \alpha)\Delta t; \ a\_7 = \alpha \Delta t. \end{aligned} \tag{38}$$

in which α, γ are defined by the assumption that the acceleration varies in each calculating step, the author selects the linear law for the varying of acceleration:

$$\ddot{\mathbf{q}}(\tau) = \ddot{\mathbf{q}}\_t + \frac{\tau}{\Delta t} (\ddot{\mathbf{q}}\_{t + \Delta t} - \ddot{\mathbf{q}}\_t) \text{with } t \le \tau \le t + \Delta t \text{ then } a = \frac{1}{2}; \gamma = \frac{1}{3}. \tag{39}$$

The condition to stabilize the roots:

$$
\Delta t \le \frac{1}{\sqrt{2}\rho\_{\text{max}}} \frac{1}{\sqrt{\alpha - \gamma}} \text{ or } \frac{\Delta t}{T\_{\text{min}}} \le \frac{1}{2\pi\sqrt{2}} \frac{1}{\sqrt{\alpha - \gamma}}\tag{40}
$$

Step 3: Calculating the stiffness matrix and the nodal force vector:

$$\mathbf{K}^\* = \overline{\mathbf{K}} + a\_0 \overline{\mathbf{M}} + a\_1 \overline{\mathbf{C}} \tag{41}$$

$$\mathbf{F}' = \overline{\mathbf{F}}\_{t + \Delta t} + \overline{\mathbf{M}} (a\_0 \mathbf{q}\_t + a\_2 \dot{\mathbf{q}}\_t + a\_3 \ddot{\mathbf{q}}\_t) + \overline{\mathbf{C}} (a\_1 \mathbf{q}\_t + a\_4 \dot{\mathbf{q}}\_t + a\_5 \ddot{\mathbf{q}}\_t) \tag{42}$$

Step 4: Determining nodal displacement vector **q***t*+Δ*t*:

$$\mathbf{K}^\*\_{t+\Delta t} \mathbf{q}\_{t+\Delta t} = \mathbf{F}^\*\_{t+\Delta t} \tag{43}$$

$$\mathbf{t} \Rightarrow \mathbf{q}\_{t + \Delta t} = \left(\mathbf{K}^\* \mathbf{r}\_{t + \Delta t}\right)^{-1} \mathbf{F}^\*\_{t + \Delta t} \tag{44}$$

repeating the loop until the time runs out.

### **3. Numerical Results of Free Vibration Analysis of Three-Layer Composite Shells with Shear Connectors**

### *3.1. Accuracy Studies*

Consider a double-curved composite shell (*0<sup>0</sup>*/*90<sup>0</sup>*/*0<sup>0</sup>*) with geometrical parameters *a* = *b*, radii *Rx* = *Ry* = *R*, thickness *h*; physical parameters *E1* = 25*E2*, *G23* = 0.2*E2*, *G13* = *G12* = 0.5*E2*, the Poisson's ratio ν12 = 0.25, and the specific weight ρ. In this case, the shear coefficient of the stub has a very large value, and this time the three-layer composite shell becomes a normal composite shell without any relative movements. We examine the convergence of the algorithm with different meshes and the comparative results of the first non-dimensional free vibration ω = ω1 *a*2*h* ρ*E*2 with Reddy [32] are shown in Table 1.


**Table 1.** The first non-dimensional fundamental frequencies with different meshes.

From Table 1 we can see clearly that, in comparison between this work and the analytical method [32], we have good agreement, demonstrating that our proposed theory and program are verified for the free vibration problem and convergence is guaranteed with 8 × 8 meshes.

### *3.2. E*ff*ects of Some Parameters on Free Vibration of the Shell*

We now consider a three-layer composite shell with geometrical parameters: length *a* is constant, width *b*, radii *Rx* = *Ry* = *R*, the total thickness *h*, the thickness of the middle layer *hc*, the thicknesses of the other layers *ht* = *hb* (*h1* = *h2* = *ht*/2, *h3* = *h4* = *hc*/2, *h5* = *h6* = *hb*/2); physical parameters: the elastic modulus *Ec* = 70 GPa, *Et* = *Eb* = 200 GPa, the Poisson's ratio ν*t* = ν*c* = <sup>ν</sup>*b* = 0.3, the specific weight ρ*c* = 2300 kg/m3,ρ*t* = ρ*b* = 7800 kg/m3, the shear coefficient of the shear connector *ktc* = *kc* = *ks*, and the shell structure is fully supported. We conduct an investigation into the first non-dimensional free vibration of the shell with non-dimensional frequencies as defined by:

$$
\overline{\omega} = a \sqrt{\frac{a^2}{h\_0}} \sqrt{\frac{\rho\_t}{E\_t}} \text{with } h\_0 = \frac{a}{50} \tag{45}
$$

3.2.1. Effect of Thickness *h*

Firstly, to examine the effect of length-to-high ratio *<sup>a</sup>*/*h*, *a* is fixed, we consider three cases with *a*/*h* = 75, 60, 50, 25, 10 (respectively). In each case, the radius-to-length ratio *R*/*a* changes from 1 to 10 as

we can see in Table 2, *b* = *a*, *hc*/*ht* = 2, and the shear coe fficient of stub *ks* = 50 MPa. The results are presented in Table 2.


**Table 2.** Effect of thickness *h* on non-dimensional fundamental frequencies.

Table 2 demonstrates that when the length-to-high ratio *a*/*h* decreases, that means the sti ffness of the structure is enhanced, correspondingly with each case of the radius-to-length ratios *R*/*<sup>a</sup>*, the non-dimensional fundamental frequency increases.

3.2.2. E ffect of the *hc*/*ht* Ratio (*ht* = *hb*)

Next, in order to study the e ffect of the *hc*/*ht* ratio, we consider five cases with *hc*/*ht*, respectively given values from 2, 4, 8, 20, 30, *b* = *a* (a is fixed), the total thickness *h* = *a*/50, and the shear coe fficient of the stub *ks* = 50 MPa. The numerical results are shown in Table 3.


**Table 3.** Effect of *hc*/*ht* ratio on non-dimensional fundamental frequencies.

Table 3 gives us a discussion that when increasing *hc*/*ht* ratio, and for *h* is constant, that means the thickness of the middle layer increases, correspondingly each case of *R*/*a* ratios, the non-dimensional fundamental frequency increases. This shows that when the thickness of the shell is constant, *hc*/*ht* increases, thus, the non-dimensional fundamental frequency increases.

### 3.2.3. E ffect of the Length-to-Width Ratio *a*/*b*

In this small section, we continually evaluate the e ffect of the length-to-width ratio *a*/*b* (a is fixed), and we meditate three cases by letting *a*/*b* = 0.5, 0.75, 1, 1.75, and 2, respectively. The total thickness of the shell *h* = *a*/50, *hc*/*ht* = 2, the radius-to-length *R*/*a* also varies from 1 to 10, as we can see in Table 4, and the shear coe fficient of stub *ks* = 50 MPa. The numerical results are tabulated in Table 4.

In Table 4 we can see obviously that, with each value of radius-to-length *R*/*<sup>a</sup>*, if the length-to-width *a*/*b* increases, the non-dimensional fundamental frequency also increases, correspondingly. This interesting point demonstrates that the sti ffness of the structure is enhanced.


**Table 4.** Effect of length-to-width ratio *a*/*b* on non-dimensional fundamental frequencies.

### 3.2.4. Effect of the Shear Coefficient of Stub ks

Finally, in this section, to examine how the shear coefficient of the stub affects the non-dimensional fundamental frequencies of this structure, we consider three cases of shear coefficient as in Table 5, and *a* = *b*, *h* = *a*/50, *hc*/*ht*= 2, *Ec* = 70 GPa is fixed. The numerical results are shown in this table.


**Table 5.** Effect of shear coefficient of the stub *ks* on non-dimensional fundamental frequencies.

In Table 5 we can see clearly that, with one value of radius-to-length *R*/*<sup>a</sup>*, when the shear coefficient of stub increases, the non-dimensional fundamental frequency of the structure ge<sup>t</sup> larger. This explains that the increasing of the shear coefficient removes the slip among layers, leading to an increase of the total stiffness of the shell structure.

### **4. Numerical Results of Forced Vibration Analysis of Three-Layer Composite Shells with Shear Connectors**

### *4.1. Accuracy Studies*

Considerign that a fully-clamped square plate with parameters can be found in [33], *a* = *b* = *1m*, *h*/*a* = 10. Material properties are the elastic modulus *E* = 30 GPa, the Poisson's ratio ν = 0.3, ρ = 2800 kg/m3. The structure is subjected to distribution sudden load *p*0 = 10<sup>4</sup> Pa. The non-dimensional displacement is calculated by the formula *w*<sup>∗</sup> = 100*Eh*<sup>3</sup> <sup>12</sup>*p*0*a*<sup>4</sup>(<sup>1</sup>−ν<sup>2</sup>)*<sup>w</sup>*0. By taking the shear coefficient and radii of the shell as very large, the comparative deflection of the centroid point of the plate between our work and [33] is shown in Figure 4, where the integral time is 5 ms, and the acting time of load is 2 ms.

We can see from Figure 4 that the deflection of the centroid point of the plate is compared to [33] is similar both shape and value. This proves that our program is verified.

**Figure 4.** The deflection of the centroid point of the plate overtime.

### *4.2. E*ff*ect of Some Parameters on the Forced Vibration of the Shell*

Now, to study effects of some parameters on forced vibration of shell, we consider a three-layer composite shell with geometrical parameters: length *a* =1 m, width *b*, thickness *h*, radii of the shell *Rx* = *Ry* = *R*, the thickness of middle layer *hc*, the thickness of other layers *ht* = *hb*. Material properties are the elastic modulus *Ec* = 8 GPa, *Et* = *Eb* = 12 GPa, the Poisson's ratio ν*t* = ν*c* = <sup>ν</sup>*b* = 0.2, the specific weight ρ*c* = 700 kg/m3, ρ*t* = ρ*b* = 2300 kg/m3, and the shear coefficient of stub *ktc* = *kcb* = *ks*. The shell is fully clamped with the uniform load *p*(*t*) varying overtime acting perpendicularly on the shell surface.

$$p(t) = \Delta P\_{\Phi} F(t); F(t) = \begin{cases} 1 - \frac{t}{\tau\_{\rm{hd}}} (0 \le t \le \tau\_{\rm{hd}})\\ 0 \text{ otherwise} \end{cases} \text{ with} \begin{cases} \Delta P\_{\Phi} = 0.20679.10^{6} \text{ N/m}^2\\ \tau\_{\rm{hd}} = 0.028 \text{ s} \end{cases} \tag{46}$$

The non-dimensional deflection and velocity of the centroid point over time are given as follows:

$$\begin{array}{l} w^\* = \frac{100h\eta^3 E\_t}{\Delta P\_\oplus a^4} w \Big( \frac{a}{2}, \frac{b}{2} \Big); \ v^\* = \frac{T h \eta^3 E\_t}{\Delta P\_\oplus a^4} v \Big( \frac{a}{2}, \frac{b}{2} \Big) \\\ u^\*\_c = \frac{10h\_0^3 E\_c}{M g u^2 (1 - \nu\_c^2)} u \Big( \frac{a}{2}, \frac{b}{2}, -\frac{h}{2} \Big); v^\*\_c = \frac{10h\_0^3 E\_c}{M g u^2 (1 - \nu\_c^2)} v \Big( \frac{a}{2}, \frac{b}{2}, -\frac{h}{2} \Big) \text{with } h\_0 = \frac{a}{50}; \ T = 0.15(s) \end{array} \tag{47}$$

where *wa*2, *b*2and *va*2, *b*2are the deflection and velocity of the centroid point of the shell.

### 4.2.1. Effect of the Length-to-High Ratio *a*/*h*

In this first small section, we study the effect of the length-to-high ratio *<sup>a</sup>*/*h*. We consider a shell with geometrical parameters *a* = *b* (a is fixed), *hc*/*ht* = 2, *R*/*a* = 6, and *a*/*h* gets value 75, 60, 50, 40 and 25, respectively, the shear coefficient of stub *ks* = 50 MPa. The non-dimensional deflection and velocity of the centroid point of the shell are presented in Figure 5 and the maximum value is shown in Table 6.

From Figure 4 and Table 6 we can see that when reducing the value of *a*/*h* ratio, this means the thickness of the shell gets thicker, the non-dimensional deflection and velocity of the centroid point overtime decrease. This is a good agreement, the reason is when the thickness of the shell increases, the stiffness of the shell obviously becomes higher.

**Figure 5.** Effect of length-to-high ratio *a*/*h* on the non-dimensional deflection and velocity of the centroid point. (**a**) Nondimensional deflection w\* versus time; and (**b**) nondimensional velocity v\* versus time.



### 4.2.2. Effect of the *hc*/*ht* Ratio (*ht* = *hb*)

Next, to investigate the effect of *hc*/*ht* ratio, we dissect the shell with geometrical parameters *a* = *b* (a is fixed); *h* = *a*/50, the value of *hc*/*ht* ratio is given as 2, 6, 8, 10, 20, and 30, *R*/*a* = 6, and the shear coefficient of stub *ks* = 50 (MPa). The non-dimensional deflection and velocity of the centroid point of the shell are shown in Figure 6, the maximum value is listed in Table 7.

We can see in Figure 6 and Table 7 that when the *hc*/*ht* ratio increases (*h* is constant), the thickness of the middle layer increases in comparison to the other layers, and the non-dimensional deflection and velocity of the centroid point overtime decreases quickly in a range from 2–20. In a range from 20–30 the non-dimensional deflection and velocity of the centroid point overtime are almost not changed. The reason is explained that when the value of the *hc*/*ht* ratio increases, the structure can reduce the ability to oscillate, and the middle layer becomes "softer", so that it imbues the vibration better than a homogeneous shell with same geometrical and physical parameters. For this particular problem, we should select the value of *hc*/*ht* ratio in a range from 20–30.

**Table 7.** Effect of *hc*/*ht* ratio on the non-dimensional deflection and velocity of the centroid point.


**Figure 6.** Effect of *hc*/*ht* ratio on the non-dimensional deflection and velocity of the centroid point. (**a**) Nondimensional velocity w\* versus time; (**b**) Nondimensional velocity v\* versus time; (**c**) Nondimensional deflection u∗*c* versus time; (**d**) Nondimensional deflection v∗*c* versus time.

### 4.2.3. Effect of the Length-to-Width Ratio *a*/*b*

We examine the effect of length-to-width ratio *a*/*b* on the non-dimensional deflection and velocity of the centroid point of the shell with *a* is fixed, *a*/*b* gets from 0.5, 1, 1.5, 2. The geometrical parameters are *h* = *a*/50, *hc*/*ht* = 2, *R*/*a* = 6, and the shear coefficient of stub *ks* = 50 MPa. The numerical results of non-dimensional deflection and velocity of the centroid point of the shell are shown in Figure 7, and the maximum value is listed in the following Table 8.

**Table 8.** Effect of the length-to-width ratio *a*/*b* on the non-dimensional deflection and velocity of the centroid point.


Now we can see in Figure 7 and Table 8, when increasing the *a*/*b* ratio, the non-dimensional deflection and velocity of the centroid point overtime decrease. This demonstrates that the stiffness of the shell gets larger, especially when the *a*/*b* ratio equals 2. This can be understood obviously that as the shape of structure gets smaller, with the same boundary condition and other parameters, the structure will become stronger.

**Figure 7.** Effect of length-to-width ratio *a*/*b* on the non-dimensional deflection and velocity of the centroid point. (**a**) Nondimensional deflection w\* versus time; and (**b**) nondimensional velocity v\* versus time.

### 4.2.4. Effect of the Shear Coefficient of the Stub

Finally, we conduct a study on the effect of the shear coefficient of the stub on the non-dimensional deflection and velocity of the centroid point of the shell. We consider four cases with *ks* = 103, 105, 1010, 1012, and 10<sup>15</sup> Pa. Geometrical parameters are *a* = *b*; *h* = *a*/50, *hc*/*ht* = 2, *R*/*a* = 6. The numerical results of non-dimensional deflection and velocity of the centroid point of the shell are plotted in Figure 8, the maximum value is shown in Table 9.

**Figure 8.** Effect of shear coefficient of stub on the non-dimensional deflection and velocity of the centroid point. (**a**) Nondimensional deflection w\* versus time; and (**b**) nondimensional velocity v\* versus time.

**Table 9.** Effect of shear coefficient of stub on the non-dimensional deflection and velocity of the centroid point.


In this last case, we can see in Figure 8 and Table 9 that when the shear coefficient of the stub increases, the non-dimensional deflection and velocity of the centroid point of the shell is reduced. It is easily understood that the enhancing of the stiffness of the stud makes the total structure ge<sup>t</sup> stronger, meaning the stiffness of the shell is increased, correspondingly.

### 4.2.5. Influence of the Mass Density of the Core Layer

Let us consider a four-edge simply supported (SSSS) shell (*b* = *a*) with *hc* = h/2, *ht* = *hb* = h/4. The shear modulus of the shear connector is *ks* = 50 MPa. The mass densities of the three layers are ρ*t* = ρ*b* = 2300 kg/m<sup>3</sup> and ρ*c* = 700, 1000, 1500, 2000, 2300 kg/m3. Nondimensional deflection and velocity of the shell center point are shown in Figure 9, maximum deflections and velocities of the shell center point are illustrated in Table 10. The mass ratio of the shell corresponding to the different values of ρ*c* compared to case ρ*t* = ρ*c* = ρ*b* = 2300 kg/m<sup>3</sup> is shown in Table 11.

**Figure 9.** Dynamic deflections of center point of the plate versus time for different ρc. (**a**) Nondimensional deflection w\* versus time, and (**b**) nondimensional velocity v\* versus time.


**Table 10.** Maximum deflections, velocities and stress of the shell center point versus time for different ρ*<sup>c</sup>*.

**Table 11.** The mass ratio of the shell corresponding to the different values of ρ*<sup>c</sup>*.


**Comment**: From the Figure 9 and Tables 10 and 11, we obtain that when the mass density of the core-layer is increased from 700 to 2300 kg, deflection and velocity of the shell center point are almost not changed. Therefore, in order to reduce the mass of the shell, we can use the triple-layer shell with shear connectors, which the core layer has a smaller mass density than other layers. Specifically, corresponding to a difference of mass density of the core layer ρ*c* = 700, 1000, 1500, 2000 kg/m3, the mass of the shell decreases by 34.79, 28.27, 17.40, and 6.56%.

### 4.2.6. Influence of Modulus of Elasticity

Let us consider a fully simply-supported (SSSS) shell (*b* = *a*) with *hc* = h/2, *ht* = *hb* = h/4. The shear modulus of the shear connector is *ks* = 50 MPa. The modulus of elasticity of the three layers are *Et* = *Eb* = 12 GPa and *Ec* = 8, 9, 10, 12 GPa. Nondimensional deflection and velocity of the shell center point are shown in Figure 10, and the maximum deflections and velocities of the shell center point are illustrated in Table 12.

**Figure 10.** Dynamic deflections of center point of the shell over time with different *Ec.*(**a**): Nondimensional deflection w\* versus time, and (**b**) nondimensional velocity v\* versus time.


**Table 12.** Maximum deflection and velocity of the shell center point over time for different *Ec.*

Comment: From the Figure 10 and Table 12, we can find that when modulus of elasticity of the core-layer is increased in a range from 8 to 12 GPa, deflection and velocity of the shell center point are slightly decreased.
