A Finite Element Model for Dynamic Analysis of Triple-Layer Composite Plates with Layers Connected by Shear Connectors Subjected to Moving Load

Triple-layered composite plates are created by joining three composite layers using shear connectors. These layers, which are assumed to be always in contact and able to move relatively to each other during deformation, could be the same or different in geometric dimensions and material. They are applied in various engineering fields such as ship-building, aircraft wing manufacturing, etc. However, there are only a few publications regarding the calculation of this kind of plate. This paper proposes novel equations, which utilize Mindlin’s theory and finite element modelling to simulate the forced vibration of triple-layered composite plates with layers connected by shear connectors subjected to a moving load. Moreover, a Matlab computation program is introduced to verify the reliability of the proposed equations, as well as the influence of some parameters, such as boundary conditions, the rigidity of the shear connector, thickness-to-length ratio, and the moving load velocity on the dynamic response of the composite plate.


Introduction
In recent years, composite layered systems have attracted a great of interest from many researchers due to their optimized material configuration, such as their high strength-to-weight ratio and stiffness-to-weight ratio. For a multi-layer composite beam, all of its layers are connected by shear connectors, which play a very crucial role in the mechanical behavior of the beam. Newmark et al. [1] proposed the governing differential equations for elastically connected steel-concrete beams, based on the linear elastic Euler-Bernoulli beam theory. Various later works [2][3][4] improved the shear effect of Newmark's model by using Timoshenko composite beam theory. Silva et al. [5] and Nguyen et al. [6,7] utilize the finite element method (FEM) to analyze the linear static properties of multi-layered composite beams. Later on, Newmark's model was developed for dynamic and non-linear problems [8][9][10][11][12]. The plate includes three layers: top layer (t), bottom layer (b) and intermediate layer (c), which can be made of the same or different materials. These layers can slide relative to each other but are not allowed to disconnect during the deformation. Each layer is attached with the local coordinate oxyt, oxyc and oxyb. The plate is divided into six components h1, h2, h3, h4, h5, h6, as in Figure 1. ut0, uc0 and ub0 are the displacements along the x axis; vt0, vc0 and vb0 represent the displacements along the y axis at the mid-plane of each layer. utc, ucb are the relative displacements between (t) layer and (c) layer along x axis; vtc, vcb are the relative displacements between (t) layer and (c) layer along y axis

Dynamic Equation for Plate Element
According to Mindlin's plate theory, the displacements of the u, v, w of each layer are shown as:  (1) Relative displacements between layers in contact surfaces are: +Relative displacement between t and c layers At the tangential plane between layers, we have:

Dynamic Equation for Plate Element
According to Mindlin's plate theory, the displacements of the u, v, w of each layer are shown as: Relative displacements between layers in contact surfaces are: +Relative displacement between t and c layers +Relative displacement between c and b layers At the tangential plane between layers, we have: where h 4 = h 3 = h c 2 . Then, we obtain: Relation between the strains and displacements of the layers is shown as: +For the k th layer: Equation (7) is rewritten in the matrix form: where: Relation between the stress and strain of the k th layer is: where: ν k is the Poisson's ratio of the k th layer and where E k is the elasticity modulus of the k th layer. We used an 8-node isoparametric element with 13 degrees of freedom (DOF) for each node. The DOFs of the i th node q i e and the plate element {q e } are defined as: 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) where Then, we have the strain of elements in the k th layer as: where B 0 ki , B 1 ki , S k are defined as: Note that B 0 ki , B 1 ki , S ki are shown in Appendix B The elastic force of shear connectors per length unit is: +For t and c layers where: in which: +For b and c layers where: in which: In the above equations, k tc ; k cb are the shear rigidity coefficients of shear connectors per length unit. We applied the principle of virtual work to the forces applied to the plate elements: By substituting Equations (1), (15), (17) and (20) into Equation (23), we obtained the dynamic equation of the plate element as follows: M eq e + K e q e = F e (t) (24) Materials 2019, 12, 598 6 of 19 With in which where L k can be seen in Appendix C With N w j = 0 0 0 0 0 0 0 0 0 0 0 0 N j (30) In the case of taking into account structural damping, we have the force vibration equation of the plate 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 the Rayleigh drag coefficients defined in [32].

Formulation of the Nodal Element Load Vector
In general cases, one considers the moving load as: particle m moves on a plate element with a non-constant velocity v in a known trajectory; the load Q is applied on the moving particle in the direction perpendicular to the plane of the plate element, as shown in Figure 1.
Let w(x,y,t) be the bending deflection of the plate under the moving load with mass m (kg). The force applied on the moving load at position (x = ξ; y = η) is [32]: where: , which is the absolute acceleration in the z direction at the position suffering moving load, is shown via the displacement vector, velocity and acceleration of the node as: x . y ∂ 2 N ml ∂x∂y + ..
Equivalent distributed force is specified according to the Delta-Dirac function as [32]: x . y ∂ 2 N ml ∂x∂y + ..
The nodal force vector of the element is computed from the distributed force p(x,y,t) applied on the element as following: where: with N ml described in Appendix D.

Differential Equation of a Triple-Layer Plate under a Moving Load
By substituting Equation (37) into Equation (31), we obtained the dynamic equation of the plate element as follows: They are linear differential equations, which have the coefficient depending on time. In order to solve these equations, we used the Newmark-beta method [32].

Accuracy Study
Example 1: in this example, the accuracy and efficiency of this approach for multi-layer plates are confirmed by comparison with the exact three-dimensional elasticity solution. The composite plate (0 • -90 • -0 • ) was subjected to the sinusoidal loading P = P 0 sin(πx/a)sin(πy/b), the material properties are E 1 = 25E 2 , G 12 = G 13 = 0.5E 2 , G 23 = 0.2E 2 , ν 12 = ν 13 = ν 23 = 0.25. The non-dimensional displacement was calculated as w c = 100h 3 E 2 P 0 a 4 w( a 2 ; b 2 ; 0). The numerical results with different meshes were compared with analytical and semi-analytical solutions, as shown in Table 1. From these results, we have shown a verification of the proposed method with six different meshes (from 10 × 10 to 20 × 20 meshes).
with the thickness of h (the thickness of each layer is h/4). The material properties and moving load parameters are the same as in Example 1. The plate was fully simply supported. The authors used the calculation program by considering the top layer and the bottom layer as two composite layers with fiber angles of 0 • and 90 • , respectively; the core layer becomes the 2-layer composite plate (90 • -0 • ); and the shear coefficient of the shear connectors was large enough. Then, the 3-layer composite plate with shear connectors became a 4-layer composite plate. The non-dimensional displacement and the non-dimensional stress (of the center point of the structure) were computed using the following The numerical results compared with [36] (based on analytical solution of layer wise theory), are presented in Table 2. From this example, we can see that the result of this work is in agreement with the results of the layer wise theory. Example 3: consider a three-layer (0 • -90 • -0 • ) composite plate with a length of a = b; a thickness of h = a/20; a modulus of elasticity of E 1 = 144.84 GPa, E 2 = 9.65 GPa, G 12 = G 13 = 4.136 GPa, G 23 = 3.447 GPa; a Poisson's ratio of ν 12 = ν 13 = ν 23 = 0.25; and a mass density of ρ = 1389.297 kg/m 3 . The plate was subjected to a moving concentrated force of P 0 = 100 N with a constant speed of v 0 = 40 m/s along the line segment y = b/2, as in reference [34]. The plate was simply supported at all edges. We verified our established program on the triple-layer plate, each layer had the thickness of 1/3 plate thickness and the shear coefficient of the shear connectors was very high. In this case, the triple-layer plate was considered as a three-layer composite, the result of this method (using the mesh 16 × 16) was compared with reference [34], as in Figure 2.
subjected to a moving concentrated force of P0 = 100 N with a constant speed of v0 = 40 m/s along the line segment y = b/2, as in reference [34]. The plate was simply supported at all edges. We verified our established program on the triple-layer plate, each layer had the thickness of 1/3 plate thickness and the shear coefficient of the shear connectors was very high. In this case, the triple-layer plate was considered as a three-layer composite, the result of this method (using the mesh 16 × 16) was compared with reference [34], as in Figure 2. The non-dimensional center deflection of the plate is defined as: Example 4: consider a simple supported rectangular plate with the pinned-free-pinned-free (SFSF-two short edges are pinned and two long edges are free) boundary condition with a length of a = 1 m; a width of b = a/2; a thickness of h = a/100; a modulus of Elasticity of E = 206.8 GPa; a Poisson's ratio of ν = 0.29; a mass density of ρ = 7820 kg/m 3 under a moving load with a mass of M = 2.3 kg along y = b/2, as in reference [33]. We then verified our established program on the triple-layer plate with three layers having the same material parameters. Each layer had the thickness of 1/3 plate thickness and the shear coefficient of the shear connector was very large. In this case, the triple-layer plate was considered as isotropic, and the result was compared with reference [33], as in Figure 3. The non-dimensional center deflection of the plate is defined as: Example 4: consider a simple supported rectangular plate with the pinned-free-pinned-free (SFSF-two short edges are pinned and two long edges are free) boundary condition with a length of a = 1 m; a width of b = a/2; a thickness of h = a/100; a modulus of Elasticity of E = 206.8 GPa; a Poisson's ratio of ν = 0.29; a mass density of ρ = 7820 kg/m 3 under a moving load with a mass of M = 2.3 kg along y = b/2, as in reference [33]. We then verified our established program on the triple-layer plate with three layers having the same material parameters. Each layer had the thickness of 1/3 plate thickness and the shear coefficient of the shear connector was very large. In this case, the triple-layer plate was considered as isotropic, and the result was compared with reference [33], as in Figure 3. From Figures 2 and 3, we can recognize that deflections of center points in our work are similar to references [33] and [34] in both shape and magnitude. This proves the reliability of our program.

Numerical Results of the Dynamic Analysis of Triple-Layer Composite Plates with Shear Connectors
Example 3: consider a triple-layer composite rectangular plate with following parameters: fixed length a = 10 m; width b, thickness h = a/50; thickness of the middle layer hs; thickness of the other two From Figures 2 and 3, we can recognize that deflections of center points in our work are similar to references [33] and [34] in both shape and magnitude. This proves the reliability of our program.

Numerical Results of the Dynamic Analysis of Triple-Layer Composite Plates with Shear Connectors
Example 3: consider a triple-layer composite rectangular plate with following parameters: fixed length a = 10 m; width b, thickness h = a/50; thickness of the middle layer h s ; thickness of the other two layers h a = h c ; modulus of elasticity of the three layers E c = 8 GPa, E t = E b = 12 GPa; Poisson's ratios ν c = 0.33; ν t = ν b = 0.2; and mass densities of three layers ρ c = 700 kg/m 3 ; ρ t = ρ b = 2300 kg/m 3 ; shear coefficients k tc = k cb = k s of the shear connector under moving load of mass M = 2.3 kg. Then, the non-dimensional deflections along z direction of the plate center point are: with h 0 = a/50 (m) and T = 2 (s).

Influence of the Modulus of Rigidity of Shear Connectors
We tested a fully simply supported (SSSS) square plate (b = a). The moving load moved along y = b/2 with velocity v = 5 m/s. The thicknesses of the three layers was h c = h/2; h t = h b = h/4. Shear moduli of the shear connector were k s = 10 2 , 10 4 , 10 6 , 10 8 , 10 10 , 10 12 , 10 14 , 10 16 Pa. Non-dimensional deflection velocity and stress of the plate center point are shown in Figure 4, maximum deflections and velocities of the plate center point are illustrated in Table 3. Nondimensional stress of the plate center point when the load moves to there is shown in Figure 5.    From Figure 4 and Table 3, we can conclude that when the shear modulus of the shear connector is increased from 10 2 to 10 14 Pa, the deflection and velocity of the plate center point, and the stress of the center point of the structure as a function of thickness in z-direction, decreased and decreased significantly at ks = 10 10 Pa (above 10 10 Pa, the deflection was nearly constant). With ks in the interval (10 2 -10 6 ), the deflection reached its maximum. We can easily see the stress jumping at the contact surfaces in Figure 5. When ks was small, the stress jumping was high. Therefore, depending on different cases in practice, we can choose correspondingly the suitable shear modulus for reducing the vibration of the plate.

Influence of the Ratio hc/ht
In the following experiment, we tested a fully simply supported (SSSS) square plate (b = a). The shear modulus of the shear connector was ks = 50 MPa. The moving load travelled along y = b/2, with velocity v = 5 m/s, hc/ht = 1, 2, 4, 6, 8, 10, 12, 14 considered (the plate thickness h = hb + hc + ht was fixed). Non-dimensional deflection and the velocity of the plate center point are shown in Figure 6. Maximum deflections and velocities of the plate center point are illustrated in Table 4. From Figure 4 and Table 3, we can conclude that when the shear modulus of the shear connector is increased from 10 2 to 10 14 Pa, the deflection and velocity of the plate center point, and the stress of the center point of the structure as a function of thickness in z-direction, decreased and decreased significantly at k s = 10 10 Pa (above 10 10 Pa, the deflection was nearly constant). With k s in the interval (10 2 -10 6 ), the deflection reached its maximum. We can easily see the stress jumping at the contact surfaces in Figure 5. When k s was small, the stress jumping was high. Therefore, depending on different cases in practice, we can choose correspondingly the suitable shear modulus for reducing the vibration of the plate.

Influence of the Ratio h c /h t
In the following experiment, we tested a fully simply supported (SSSS) square plate (b = a).   Figure 6 and Table 4, we can see that when the ratio hc/ht was increased from 8 to 14 (i.e., increasing the thickness of the middle layer since the plate thickness h is fixed), deflection and velocity of the plate center point decreased significantly. This proved that the middle layer, which had smaller elastic modulus than the other two layers, could absorb vibration better and thus reduce deflection and velocity of the plate center point as a function of thickness in z-direction. Therefore, we proposed that we do not have to use material with large modulus of elasticity to decrease the vibration. Instead, we can use a triple-layer plate (with a shear connector) in which the middle layer has a smaller elastic modulus with a specific thickness ratio. Here, we obtain an interesting property: the vertical displacement w decreased when the ratio hc/ht increased, and in contrast, the horizontal displacements uc and vc (at the contact surfaces) increased. This can be explained by the strain energy focusing on the membrane direction.  From Figure 6 and Table 4, we can see that when the ratio h c /h t was increased from 8 to 14 (i.e., increasing the thickness of the middle layer since the plate thickness h is fixed), deflection and velocity of the plate center point decreased significantly. This proved that the middle layer, which had smaller elastic modulus than the other two layers, could absorb vibration better and thus reduce deflection and velocity of the plate center point as a function of thickness in z-direction. Therefore, we proposed that we do not have to use material with large modulus of elasticity to decrease the vibration. Instead, we can use a triple-layer plate (with a shear connector) in which the middle layer has a smaller elastic modulus with a specific thickness ratio. Here, we obtain an interesting property: the vertical displacement w decreased when the ratio h c /h t increased, and in contrast, the horizontal displacements u c and v c (at the contact surfaces) increased. This can be explained by the strain energy focusing on the membrane direction.

Influence of the Ratio a/h
Let us consider a fully simply supported (SSSS) square plate (b = a) with h c = h/2, h t = h b = h/4. Shear modulus of the shear connector was k s = 50 MPa. The moving load travelled along y = b/2 with velocity v = 5 m/s. a/h = 30, 40, 50, 60, 70, 80, 90, 100 used. Non-dimensional deflection and the velocity of the plate center point are shown in Figure 7, maximum deflections and velocities of the plate center point are illustrated are Table 5.   Figure 7 and Table 5, we can recognize that when the plate thickness was increased from a/100 to a/30, deflection, velocity and stress of the plate center point decreased and decreased significantly in the range of a/50 to a/30.    Figure 7 and Table 5, we can recognize that when the plate thickness was increased from a/100 to a/30, deflection, velocity and stress of the plate center point decreased and decreased significantly in the range of a/50 to a/30. From Figure 8 and Table 6, we can find that when the velocity of the moving load was increased from 5 to 40 m/s, deflection, velocity of the plate center point as a function of thickness in z-direction  Table 7.   From Figure 8 and Table 6, we can find that when the velocity of the moving load was increased from 5 to 40 m/s, deflection, velocity of the plate center point as a function of thickness in z-direction  Table 7.  Table 7.   From Figure 9 and Table 7, we can find that when the mass density of the core-layer was increased from 700 to 2300 kg, deflection and velocities of the plate center point as a function of thickness in z-direction were almost not changed. Therefore, in order to reduce the mass of the plate we can use a triple-layer plate with shear connectors, where the core layer has a smaller mass density than other layers.

Influence of Modulus of Elasticity
Let us consider a fully simply supported (SSSS) square plate (b = a) with h c = h/2, h t = h b = h/4. The shear modulus of the shear connector was k s = 50 MPa. The moving load travelled along y = b/2 with velocity v = 5m/s. A modulus of elasticity of three layers E t = E b = 12 GPa and E c = 8, 9, 10, 12 GPa, was used. Non-dimensional deflection and velocity of the plate center point is shown in Figure 10, maximum deflections and velocities of the plate center point are illustrated in Table 8.   Figure 9 and Table 7, we can find that when the mass density of the core-layer was increased from 700 to 2300 kg, deflection and velocities of the plate center point as a function of thickness in z-direction were almost not changed. Therefore, in order to reduce the mass of the plate we can use a triple-layer plate with shear connectors, where the core layer has a smaller mass density than other layers.

Influence of Modulus of Elasticity
Let us consider a fully simply supported (SSSS) square plate (b = a) with hc = h/2, ht = hb = h/4. The shear modulus of the shear connector was ks = 50 MPa. The moving load travelled along y = b/2 with velocity v = 5m/s. A modulus of elasticity of three layers Et = Eb = 12 GPa and Ec = 8, 9, 10, 12 GPa, was used. Non-dimensional deflection and velocity of the plate center point is shown in Figure 10, maximum deflections and velocities of the plate center point are illustrated in Table 8.   Figure 10 and Table 8, we can find that when modulus of elasticity core-layer was increased from 8 to 12 GPa, deflection, velocity of the plate center point and the stress of center point as a function of thickness in z-direction were slightly decreased.

Conclusions
By using published theories of multi-layered beams, Mindlin's plate theory and finite element modelling, we simulated the forced vibration of a triple-layer plate with layers connected by shear connectors subjected to a moving load. The influences of some structural parameters on the dynamic response of the plate were also examined in the paper. We found that, in the regarded plate, the shear  From Figure 10 and Table 8, we can find that when modulus of elasticity core-layer was increased from 8 to 12 GPa, deflection, velocity of the plate center point and the stress of center point as a function of thickness in z-direction were slightly decreased.

Conclusions
By using published theories of multi-layered beams, Mindlin's plate theory and finite element modelling, we simulated the forced vibration of a triple-layer plate with layers connected by shear connectors subjected to a moving load. The influences of some structural parameters on the dynamic response of the plate were also examined in the paper. We found that, in the regarded plate, the shear coefficient of the shear connector played a crucial role. Specifically, when the shear coefficient was sufficiently large, the plate could be considered a sandwich plate. The flexibility of the shear coefficient helped engineer a plate with the desired mechanical properties. From the results of numerical tests, we proposed that to reduce plate vibration, the elastic modulus of the middle layer should be smaller than the outside two layers, and the thickness of the middle layers should be 20-30 times larger than the two outside layers.

Conflicts of Interest:
The authors declare no conflict of interest.