Finite Element Modelling of a Composite Shell with Shear Connectors

A three-layer composite shell with shear connectors is made of three shell layers with one another connected by stubs at the contact surfaces. These layers can have similar or different geometrical and physical properties with the assumption that they always contact and have relative movement in the working process. Due to these characteristics, they are used widely in many engineering applications, such as ship manufacturing and production, aerospace technologies, transportation, and so on. However, there are not many studies on these types of structures. This paper is based on the first-order shear deformation Mindlin plate theory and finite element method (FEM) to establish the oscillator equations of the shell structure under dynamic load. The authors construct the calculation program in the MATLAB environment and verify the accuracy of the established program. Based on this approach, we study the effects of some of the geometrical and physical parameters on the dynamic responses of the shell.


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 great 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 Symmetry 2019, 11, 527 2 of 20 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 effect due to the shear locking coefficient 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 different 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 efficient 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 effect 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.

Equation of Motion of the Shell Element
Consider a three-layer composite shell with shear connectors as shown in Figure 1.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 h1, h2, h3, h4, h5, h6 as shown in Figure 1; ut0, uc0, and ub0 represent displacements in x direction; vb0 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: ( ) 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 And for layer c and layer b we have: 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; u t0 , u c0 , and u b0 represent displacements in x direction; v b0 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 (x k , y k , z k ) of layer k are expressed as follows: where ϕ k and ψ k are the transverse normal rotations of the x k and y k directions.The relative movements among the contact surfaces are defined by the following equations For the layer t and layer c we have And for layer c and layer b we have: Symmetry 2019, 11, 527 4 of 20 Note that at the contact surfaces, we have: with h 4 = h 3 = h c 2 .From Equations ( 1)-( 4), we get: The relation between strain and displacement of each layer is expressed as follows For the layer k, we have: We can rewrite in a matrix form as follow in which The relation between stress and strain of layer k is expressed as followIs necessary bild?
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.
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 i e and the total degree of freedom of the shell element q e is defined as follow.
q e = q 1 e q 2 e q 3 e q 4 e q 5 e q 6 e q 7 e q 8 e T ( 13) in which N i (i = 1 ÷ 8) can be defined as in [28].
The relation between stress and strain of layer k is expressed as followIs necessary bild?
, k k D G are the bending rigidity and shear rigidity of layer k, respectively, and k ν is the Poisson ratio of layer k.
In this work, the thickness of the shell is thin or medium ( 100 10 a a h = ÷ , 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 i e q and the total degree of freedom of the shell element e q is defined as follow.By substituting in the expression for verifying displacement of element we have: in which B 0 k ; B 1 k ; S k are defined as follows where B 0 ki , B 1 ki and S 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: With Symmetry 2019, 11, 527 6 of 20 in which For layer c and b we have with in which Here, k tc and k cb 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: By substituting Equations ( 1), ( 15), (17), and (20) into Equation ( 23), we obtain the dynamic equation of the shell element as follows: with in which where in which with N w j = 0 0 0 0 0 0 0 00000 N j (30) In the case of taking into account the structural damping, we have the force vibration equation of the shell element as follows: M e .. q e + C e .q e + K e q e = F e (t) (31) in which C e = αM e + βK e and α, β are Rayleigh drag coefficients defined in [30,31].

The Differential 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: 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.

Accuracy Studies
Consider a double-curved composite shell (0 0 /90 0 /0 0 ) with geometrical parameters a = b, radii Rx = Ry = R, thickness h; physical parameters E1 = 25E2, G23 = 0.2E2, G13 = G12 = 0.5E2, 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, For the free vibration analysis, the natural frequencies can be obtained by solving the equation: or in another form: where ω is the natural frequency.Flowchart of Newmark-beta method [31] Step 1: Determine the first conditions: .
From the first conditions, we obtain: ..
q t+∆t by q t+∆t , we have where: 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: .. q(τ) = ..
The condition to stabilize the roots: Step 3: Calculating the stiffness matrix and the nodal force vector: Step 4: Determining nodal displacement vector q t+∆t : repeating the loop until the time runs out.

Accuracy Studies
Consider a double-curved composite shell (0 0 /90 0 /0 0 ) with geometrical parameters a = b, radii 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.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.

Effects 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 R x = R y = R, the total thickness h, the thickness of the middle layer h c , the thicknesses of the other layers 3 , the shear coefficient of the shear connector k tc = k c = k s , 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:

Effect of Thickness h
Firstly, to examine the effect of length-to-high ratio a/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, h c /h t = 2, and the shear coefficient of stub k s = 50 MPa.The results are presented in Table 2. Table 2 demonstrates that when the length-to-high ratio a/h decreases, that means the stiffness of the structure is enhanced, correspondingly with each case of the radius-to-length ratios R/a, the non-dimensional fundamental frequency increases.

Effect of the h
Next, in order to study the effect of the h c /h t ratio, we consider five cases with h c /h t , respectively given values from 2, 4, 8, 20, 30, b = a (a is fixed), the total thickness h = a/50, and the shear coefficient of the stub k s = 50 MPa.The numerical results are shown in Table 3. Table 3 gives us a discussion that when increasing h c /h t 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, h c /h t increases, thus, the non-dimensional fundamental frequency increases.

Effect of the Length-to-Width Ratio a/b
In this small section, we continually evaluate the effect 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, h c /h t = 2, the radius-to-length R/a also varies from 1 to 10, as we can see in Table 4, and the shear coefficient of stub k s = 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/a, if the length-to-width a/b increases, the non-dimensional fundamental frequency also increases, correspondingly.This interesting point demonstrates that the stiffness of the structure is enhanced.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, h c /h t = 2, E c = 70 GPa is fixed.The numerical results are shown in this table.In Table 5 we can see clearly that, with one value of radius-to-length R/a, when the shear coefficient of stub increases, the non-dimensional fundamental frequency of the structure get 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.

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/m 3 .The structure is subjected to distribution sudden load p 0 = 10 4 Pa.The non-dimensional displacement is calculated by the formula w * = 100Eh 3 12p 0 a 4 (1−ν 2 ) w 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.
displacement is calculated by the formula .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.

Effect 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 0.2

Effect 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 R x = R y = R, the thickness of middle layer h c , the thickness of other layers h t = h b .Material properties are the elastic modulus E c = 8 GPa, E t = E b = 12 GPa, the Poisson's ratio ν t = ν c = ν b = 0.2, the specific weight ρ c = 700 kg/m 3 , ρ t = ρ b = 2300 kg/m 3 , and the shear coefficient of stub k tc = k cb = k s .The shell is fully clamped with the uniform load p(t) varying overtime acting perpendicularly on the shell surface.
The non-dimensional deflection and velocity of the centroid point over time are given as follows: where w a 2 , b 2 and v a 2 , b 2 are the deflection and velocity of the centroid point of the shell.

Effect of the Length-to-High Ratio a/h
In this first small section, we study the effect of the length-to-high ratio a/h.We consider a shell with geometrical parameters a = b (a is fixed), h c /h t = 2, R/a = 6, and a/h gets value 75, 60, 50, 40 and 25, respectively, the shear coefficient of stub k s = 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.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.Next, to investigate the effect of h c /h t ratio, we dissect the shell with geometrical parameters a = b (a is fixed); h = a/50, the value of h c /h t ratio is given as 2, 6, 8, 10, 20, and 30, R/a = 6, and the shear coefficient of stub k s = 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 h c /h t 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 h c /h t 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 h c /h t ratio in a range from 20-30.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 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, h c /h t = 2, R/a = 6, and the shear coefficient of stub k s = 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.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.

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 = 10 3 , 10 5 , 10 10 , 10 12 , and 10 15 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.
Table 9.Effect of shear coefficient of stub on the non-dimensional deflection and velocity of the centroid point.

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 k s = 10 3 , 10 5 , 10 10 , 10 12 , and 10 15 Pa.Geometrical parameters are a = b; h = a/50, h c /h t = 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.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 get stronger, meaning the stiffness of the shell is increased, correspondingly.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 get stronger, meaning the stiffness of the shell is increased, correspondingly.Comment: From the Figure 9 and Tables 10, 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, the core layer has a smaller mass density than other layers.Specifically, corresponding to a difference of mass density of the core layer  12.The mass ratio (100 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/m 3 , the mass of the shell decreases by 34.79, 28.27, 17.40, and 6.56%.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.

Conclusions
finite element method (FEM) is the numerical method for solving problems of engineering and mathematical physics, including the calculation of shell structures.Establishing the balance equation describing the vibration of shell structure is quite simple and it is very convenient for coding on a personal computer (PC).The proposed program is able to analyze the static bending, dynamic response, nonlinear problems, etc., with complicated structures, which are not easy to solve by analytical methods.
Based on the finite element method, we established the equilibrium equation of a triple-layer composite shell with shear connectors subjected to dynamic loads.In this paper, employing of the eight-node isometric element is suitable.To exactly describe the strain field, the displacements of the three-layer shell with shear connectors, and the 13-degrees of freedom element is used, in which the three layers have the same a degree of freedom in the z-direction, and the other 12 degrees of freedom are described as the linear displacement and rotation angle of each layer.Hence, the displacement field and the strain field of each layer can be investigated deeply.We have created the program in the MATLAB environment to investigate effects of various geometrical parameters on free and forced vibrations of shells.To sum up, some main interesting points of this paper are listed in the following statements.
In general, the geometrical parameters effect strongly on free and forced vibrations of the shell; when the shape of the shell is small, the structure gets stiffer.
Based on the numerical results, we realized that for this type of structure, the shear coefficient of the stub plays a very important role.Especially, when the stiffness of the shear coefficient is large enough, this structure seems to be a sandwich shell.
From the above computed results, we suggest that in order to reduce the vibration of such a structure, we should use the middle layer, having the elastic modulus less than other layers, and the thickness of the middle layer 20-30 times larger than the other ones.
We suggest that, in order to reduce the volume of the shell structure subjected to the blast load, we should consider the triple-layer shell with the core layer having a smaller density than the two 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.

Conclusions
The finite element method (FEM) is the numerical method for solving problems of engineering and mathematical physics, including the calculation of shell structures.Establishing the balance equation describing the vibration of shell structure is quite simple and it is very convenient for coding on a personal computer (PC).The proposed program is able to analyze the static bending, dynamic response, nonlinear problems, etc., with complicated structures, which are not easy to solve by analytical methods.
Based on the finite element method, we established the equilibrium equation of a triple-layer composite shell with shear connectors subjected to dynamic loads.In this paper, employing of the eight-node isometric element is suitable.To exactly describe the strain field, the displacements of the three-layer shell with shear connectors, and the 13-degrees of freedom element is used, in which the three layers have the same a degree of freedom in the z-direction, and the other 12 degrees of freedom are described as the linear displacement and rotation angle of each layer.Hence, the displacement field and the strain field of each layer can be investigated deeply.We have created the program in the MATLAB environment to investigate effects of various geometrical parameters on free and forced vibrations of shells.To sum up, some main interesting points of this paper are listed in the following statements.0 0 0 0 0 0 0000 N i R y ∂N i ∂y ∂N i ∂x 0 0 0 0 0 0 0000 0 0 0 0 0 ∂N i ∂y 0 0 0000 0 0 0 00000 0 N i 0 00000 0 0 0 10z t 0000000000 010z t 000000000 0000000000001 000010z c 000000 0000010z c 00000 0000000000001

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.

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.

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

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

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

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

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

3 ,
and the shear coefficient of stub ktc = kcb = ks.The

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

Figure 5 .
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.

Figure 5 .
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.

Figure 6 .
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.

Figure 6 .
Figure 6.Effect of h c /h t 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

Figure 7 .
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.

Figure 7 .
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.

Figure 8 .
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.

4. 2 . 5 .Figure 8 .
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.

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 h c = h/2, h t = h b = h/4.The shear modulus of the shear connector is k s = 50 MPa.The mass densities of the three layers are ρ t = ρ b = 2300 kg/m 3 and ρ c = 700, 1000, 1500, 2000, 2300 kg/m 3 .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 3 is shown in

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

kg/m 3 , 4 . 2 . 6 .
the mass of the shell decreases by 34.79, 28.27, 17.40, and 6.56%.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

Figure 9 .Table 10 .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 ρ c .Maximum Values ρ c = 700 (kg/m 3 )

4. 2 . 6 .
Influence of Modulus of ElasticityLet us consider fully simply-supported (SSSS) shell (b = a) with h c = h/2, h t = h b = h/4.The shear modulus of the shear connector is k s = 50 MPa.The modulus of elasticity of the three layers are E t = E b = 12 GPa and E c = 8, 9, 10, 12 GPa.Nondimensional deflection and velocity of the shell center point are shown in Figure10, and the maximum deflections and velocities of the shell center point are illustrated in Table12 .

Figure 10 .
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.

Figure 10 .Table 12 .
Figure 10.Dynamic deflections of center point of the shell over time with different E c .(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 E c .Maximum Values E c = 8 GPa E c = 9 GPa E c = 10 GPa E c = 12 GPa E c = 12 GPa w * max 3.9020 3.7494 3.5895 3.4252 3.3031 v * max

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

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

Table 3 .
Effect of h c /h t ratio on non-dimensional fundamental frequencies.

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

Table 5 .
Effect of shear coefficient of the stub k s on non-dimensional fundamental frequencies.

Table 6 .
Effect of length-to-high ratio a/h the non-dimensional deflection and velocity of the centroid point.
4.2.2.Effect of the hc /h t Ratio (h t = h b )

Table 7 .
Effect of h c /h t ratio on the non-dimensional deflection and velocity of the centroid point.

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

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