MESHLESS MODELLING OF LAMINATE MINDLIN PLATES UNDER DYNAMIC LOADS MESHLESS MODELLING OF LAMINATE MINDLIN PLATES UNDER DYNAMIC LOADS

Collocation method and Galerkin method have been dominant in the existing meshless methods. A meshless local Petrov-Galerkin (MLPG) method is applied to solve laminate plate problems described by the Reissner-Mindlin theory for transient dynamic loads. The Reiss-ner-Mindlin theory reduces the original three-dimensional (3-D) thick plate problem to a two-dimensional (2-D) problem. The bending moment and the shear force expressions are obtained by integration through the laminated plate for the considered constitutive equations in each lamina. The weak-form on small subdomains with a Heaviside step function as the test functions is applied to derive local integral equations. After performing the spatial MLS approximation, a system of ordinary differential equations of the second order for certain nodal unknowns is obtained. The derived ordinary differential equations are solved by the Houbolt finite-difference scheme as a time-stepping method.


Introduction
Composite materials are now common engineering materials used in a wide range of applications. They play an important role in the aviation, aerospace and automotive industry, are also used in the construction of ships, submarines, nuclear and chemical facilities, etc. During the last several decades, laminated composite plates have been widely used in engineering structures. The optimization of the topology of structural lay out and composite plates and shells has great impact on the performance of structures [1][2][3].
Previous research results show that the transverse shear effects are more significant for orthotropic plates than for isotropic ones [4,5]. It is well known that the classical thin plate theory of Kirchhoff gives rise to certain non-physical simplifications mainly related to the omission of the shear deformation and the rotary inertia, which become more significant for increasing thickness of the plate. The effects of shear deformation and rotary inertia are taken into account in the Reissner-Mindlin plate bending theory [6]. A higherorder theory for plates deformed in shear and normal mode was used by Qian et al. [7]. There, the plate material is made of two isotropic constituents and exhibits macroscopically isotropic material properties which vary in the thickness direction only. Results of the static analysis and computed natural frequencies of a simply supported square plate match well with corresponding analytical values.
In addition, the governing equations for thick orthotropic shells are also quite complicated. A review on early applications of BEM to shells is given by Beskos [8]. The analysis of thin elastic plates by the boundary element method (BEM) is a research subject since many years. Much previous research works have been done for static and dynamic analysis of isotropic thin plates by the BEM [9]. The first application of the boundary integral equation method to Reissner's plate model was given by Van der Ween [10]. Dynamic analysis of elastic Reissner-Mindlin plates was performed by the direct BEM in the frequency domain [11,12]. All previous BEM applications deal with isotropic Reissner-Mindlin plates. Wang and Huang [13] were the first who applied BEM to orthotropic thick plates.
In spite of the great success of the FEM and the BEM as accurate and effective numerical tools for the solution of boundary or initial-boundary value problems in domains with complex shapes, there is still a growing interest in developing new advanced numerical methods [14,15]. Meshless approaches for problems of continuum mechanics have attracted much attention during the past decade especially owing to their high adaptivity and low costs to prepare input data for numerical analyses. Many meshless methods are derived from a weak-form formulation on global domain or a set of local subdomains [16,17]. In the global formulation background cells are required for the integration of the weak-form. In methods based on local weak-form formulation no cells are required and therefore they are often referred to as truly meshless methods. If a simple form is chosen for the geometry of the subdomains, numerical integrations over them can be easily carried out. The first application of a meshless method to plate/shell problems was given by Krysl and Belytschko [18,19], where they applied the elementfree Galerkin method. The Moving Least-Square (MLS) approximation yields a C 1 -continuity which satisfies the Kirchhoff hypotheses. The continuity of the MLS approximation is given by the minimum between the continuity of the basis functions and that of the weight function. So continuity can be tuned to a desired value. Their results showed an excellent convergence, however, their formulation is not applicable to shear deformable plate/shell problems. Recently, Noguchi et al. [20] used a mapping technique to transform a curved surface into a flat two-dimensional space. Then, the element-free Galerkin method can be applied also to thick plates or shells including the shear deformation effects. The meshless local Petrov-Galerkin (MLPG) method is a fundamental base for the derivation of many meshless formulations, since the trial and the test functions can be chosen from different functional spaces. The method has been successfully applied also to plate and shell problems with homogeneous material properties [21][22][23][24].

Theory Backround of Governing Equations
The classical laminate plate theory is an extension of the classical plate theory to composite laminates. Consider a plate of total thickness h composed of N orthotropic layers with the mean surface occupying the domain Ω in the plane (x 1 , x 2 ). The x 3 ϵ z axis is perpendicular to the mid-plane (Fig.1). The k-th layer is located between the points z ϭ z k and z ϭ z kϩ1 in the thickness direction.
The Reissner-Mindlin plate bending theory is used to describe the plate deformation. The transverse shear strains are represented as constant throughout the plate thickness and some correction coefficients are required for the computation of transverse shear forces in that theory. Then, the spatial displacement field in time τ, due to transverse loading and expressed in terms of displacement components u 1 , u 2 and u 3 , has the following form [25] where x ϭ [x 1 , x 2 ] T is the position vector, and w α (x 1 , x 2 , τ) represent the rotations around the in-plane axes and the out-of-plane deflection, respectively (Fig. 1). The linear strains are given by In the case of orthotropic materials for the k-th lamina, the relation between the stresses σ ij and the strains ε ij is described by the constitutive equations for the stress tensor where the material stiffness coefficients c (k) ijml are assumed to be homogeneous for the k-th lamina.
It can be seen from equation (2) that the strains are continuous throughout the plate thickness. Hence, discontinuous material coefficients yield discontinuities in stresses on the lamina surfaces.
For plane problems the constitutive equation (3) is frequently written in terms of the second-order tensor of elastic constants. The constitutive equation for orthotropic materials and plane stress problem are given for example in [23,26].
Despite the stress discontinuities, one can define the integral quantities such as the bending moments M αβ and the shear forces Q αβ as (4) and , where κ ϭ 5/6 in the Reissner plate theory.
Substituting constitutive equations and (2) into moment and force resultants (5) allows the expression of the bending moments M αβ and shear forces Q α for α, β ϭ 1, 2, in terms of rotations and lateral displacements of the orthotropic plate. In the case of the considered layer-wise continuous material properties through the plate thickness, one obtains In eq. (6), repeated indices α and β do not imply summation, and the material parameters D αβ and C αβ are given as (8) For a homogeneous plate equations (8) are reduced into simple forms , , , The plate is subjected to a transverse dynamic load q(x, t). Assuming the mass density to be homogeneous within each lamina and using the Reissner's linear theory of thick plates [25], the equations of motion may be written as where are global inertial characteristics of the laminate plate. If the mass density is constant throughout the plate thickness, we obtain .
Throughout the analysis, Greek indices vary from 1 to 2, and the dots over a quantity indicate differentiations with respect to time τ.

Local Petrov-Galerkin Weak-Form
Instead of writing the global weak-form for the above governing equations, the MLPG methods construct the weak-form over local subdomains such as Ω s , which is a small region taken for each node inside the global domain [17]. The local subdomains overlap each other and cover the whole global domain Ω (Fig. 2). The local subdomains could be of any geometrical shape and size. In the current paper, the local subdomains are taken to be of circular shape.

Fig. 2. Local boundaries for weak formulation, the domain Ω x for MLS approximation of the trial function, and support area of weight function around node x i
The local weak-form of the governing equations (9) and (10) for x i ʦ Ω i s can be written as (12) here w * αβ (x) and w * (x) are weight or test functions.
Applying the Gauss divergence theorem to Eqs. (11) and (12) one obtains (13)   (14) where ѨΩ i s is the boundary of the local subdomain and M α (x, τ) ϭ ϭ M αβ (x, τ)n β (x) is the normal bending moment and n α is the unit outward normal vector to the boundary ѨΩ i s . The local weak-forms (13) and (14) are the starting point for deriving local boundary integral equations on the basis of appropriate test functions. Unit step functions are chosen for the test functions w * αβ (x) and w * (x) in each subdomain , (15) .
Then, the local weak-forms (13) and (14) are transformed into the following local integral equations (LIEs) (16) In the above local integral equations, the trial functions w α (x, τ) related to rotations, and w 3 (x, τ) related to transversal displacements, are chosen as the moving least-squares (MLS) approximations over a number of nodes randomly spread within the domain of influence.

Numerical Implementation
In general, a meshless method uses a local interpolation to represent the trial function with the values (or the fictitious values) of the unknown variable at some randomly located nodes. The , , , .
, , # moving least-squares (MLS) approximation [27] used in the present analysis may be considered as one of such schemes. Let us consider a sub-domain Ω x of the problem domain Ω in the neighbourhood of a point x for the definition of the MLS approximation of the trial function around x (Fig. 3). To approximate the distribution of the generalized displacements (rotations and deflection) in Ω x over a number of randomly located nodes {x a }, a ϭ 1, 2, … n, where ] is a complete monomial basis of order m, and a ෂ (x, τ) ϭ [a 1 (x, τ), where ν a (x) Ͼ 0 is the weight function associated with the node a and the square power is considered in the sense of scalar product. Recall that n is the number of nodes in Ω x for which the weight function ν a (x) Ͼ 0 and ŵ a (τ) are the fictitious nodal values, but not the nodal values of the unknown trial function w h (x, τ), in general. The stationarity of J in eq. (19) with respect to a ෂ (x, τ) leads to The solution of eq. (20) for and the subsequent substitution into eq. (18) lead to the following expression (22) where In eq. (22), φ a (x) is usually referred to as the shape function of the MLS approximation corresponding to the nodal point x a . From eqs. (21) and (23), it can be seen that φ a (x) ϭ 0 when v a (x) ϭ 0. In practical applications, v a (x) is often chosen in such a way that it is non-zero over the support of the nodal point x a . The support of the nodal point x a is usually taken to be a circle of the radius r a centred at x a (see Fig. 3). The radius r a is an important parameter of the MLS approximation because it determines the range of the / interaction (coupling) between the degrees of freedom defined at considered nodes.
Usually a 4 th -order spline-type weight function is applied [17] where d a ϭ ʈx Ϫ x a ʈ and r a is the radius of the circular support domain. With eq. (24), the C 1 -continuity of the weight function is ensured over the entire domain, therefore the continuity condition of the bending moments and the shear forces is satisfied. The size of the support r a should be large enough to cover a sufficient number of nodes in the domain of definition to ensure the regularity of the matrix A. The value of n is determined by the number of nodes lying in the support domain with radius r a .
The partial derivatives of the MLS shape functions are obtained as [28] (25) ,k ϭ (A Ϫ1 ) ,k represents the derivative of the inverse of A with respect to x k , which is given by The directional derivatives of w(x, τ) are approximated in terms of the same nodal values as (26) Substituting the approximation (26) into the definition of the bending moments (6) and then using M α (x, τ) ϭ M αβ (x, τ)n β (x), one obtains for M(x, τ) ϭ [M 1 (x, τ), M 2 (x, τ)] T (27) where the vector w *a (τ) is defined as a column vector w *a (τ) ϭ ϭ [ŵ a 1 (τ), ŵ a 2 (τ)] T , the matrices N α (x) are related to the normal vector n(x) on ѨΩ s by and and the matrices B a α are represented by the gradients of the shape functions as , .
The influence of the material properties for composite laminates is incorporated into C αβ and D αβ defined in equations (7). Similarly one can obtain the approximation for the shear forces  Then, insertion of the MLS-discretized moment and force fields (27) and (28) into the local integral equations (16) and (17) yields the discretized local integral equations It should be noted here that there are neither Lagrange multipliers nor penalty parameters introduced into the local weakforms (12) because the essential boundary conditions on Γ i sw (part of the global boundary with prescribed rotations or displacements) can be imposed directly, using the interpolation approximation (22) for , where w ෂ (x i , τ) is the generalized displacement vector prescribed on the boundary Γ i sw . For a clamped plate all three vector components (rotations and deflection) are vanishing on the fixed edge, and eq. (31) is used at all the boundary nodes in such a case. However, for a simply supported plate only the third component of the displacement vector (deflection) is prescribed, while the rotations are unknown. Then, the entire equation (29) and the third component of eq. (31) are applied to the nodes lying on the global boundary. On those parts of the global boundary where no displacement boundary conditions are prescribed both local integral equations (29) and (30) are applied.
* a a a a a n 3 1 Collecting the discretized local boundary-domain integral equations together with the discretized boundary conditions for the generalized displacements, one obtains a complete system of ordinary differential equations and it can be rearranged in such a way that all known quantities are on the r.h.s. Thus, in matrix form the system becomes where M is the mass matrix, K is the stiffness matrix and F is the external load vector. Recall that the system matrix has a block structure. There are many time integration procedures for the solution of this system of ordinary differential equations. In the present work, the Houbolt method is applied [29]. In the Houbolt finitedifference scheme the value of the time-step has to be appropriately selected with respect to material parameters (elastic wave velocities) and time dependence of the boundary conditions.

Numerical Examples
In this section, numerical results are presented for two laminate plates under an impact load with Heaviside time variation. Clamped and simply supported square plates are analysed. In order to test the accuracy, the numerical results obtained by the present method are compared with the results provided by the FEM-ANSYS code using a very fine mesh. In numerical calculations, deviations of the results for plates with a laminate structure from those corresponding to a homogeneous plate are investigated. The mass density is assumed to be uniform m within the whole bulk of the plate.

Example 1: Clamped square plate
In this example the clamped orthotropic thick square plate under an impact load with Heaviside time variation is analyzed. The clamped square plate has side-length a ϭ 0.254m and the plate thicknesses is h/a ϭ 0.05. Homogeneous material properties are considered firstly to test the accuracy of the present computational method. The following material parameters are used in our numerical analysis: Young's moduli E 2 ϭ 0.6895 и 10 10 N/m 2 and E 1 ϭ 2E 2 , Poisson's ratios ν 21 ϭ 0.15 and ν 12 ϭ 0.3 and mass density ρ ϭ 5.0 ϫ 10 3 kgm Ϫ3 . The used shear moduli correspond to Young's modulus E 2 , namely, For the numerical modelling we used again 441 nodes with a regular distribution (Fig. 3). Numerical calculations are carried out for a time-step Δτ ϭ 0.357 и 10 Ϫ4 s. The MLPG results are compared with those obtained by FEM-ANSYS computer code, where 900 quadrilateral eight-node shell elements with 1000 time increments were used. The static central deflection is w 3 stat (0) ϭ ϭ 8.842 и 10 Ϫ3 m for the considered load q 0 ϭ 2.07 ϫ 10 6 Nm Ϫ2 . The static bending moment is M 11 stat (0) ϭ 3064 Nm. The time is normalized by τ 0 ϭ a 2 /4Ί ρh / D ϭ 0.3574 и 10 Ϫ2 s, where D ϭ ϭ Eh 3 /[12(1 Ϫ ν 2 )] is bending rigidity of plate. A good agreement of the present results for the deflection and the bending moment at the plate center and the FEM results is observed.
The peaks of the moment amplitudes are shifted to shorter time instants for the orthotropic plate with a larger flexural rigidity. Since the mass density is the same for both isotropic and orthotropic plates, the wave velocity is higher for the orthotropic plate with higher Young's modulus. The amplification of the bending moments due to the dynamic impact load for both isotropic and orthotropic plates are almost the same if they are normalized with respect to the corresponding static values. The static bending moment for the orthotropic plate is slightly higher at the center of the plate [30]. Time variation of the deflection at the center of a clamped square plate subjected to a suddenly applied load Fig.  4. Then, the peak value for the orthotropic plate under an impact load has to be higher since that value is normalized by corresponding static value at the center of the isotropic plate. A simply supported three-ply orthotropic square plate under an impact load with Heaviside time variation is analyzed. The used geometrical and material parameters are the same as for the previous clamped plate. The total plate thickness is h ϭ 0.0127m, which is equal to the thickness of the homogeneous plate in previous examples. The bottom and top layers have the same thickness h 1 ϭ h 3 ϭ h/4. Young's moduli for both bottom and top layers are the same and they are 5 times larger than the ones corresponding to the homogeneous orthotropic plate. The second mid-layer has the thickness h 2 ϭ h/2 and the same material properties as the homogeneous plate analyzed in previous example. For the numerical modelling we used again 441 nodes with a regular distribution. A uniformly distributed load q 0 ϭ 300 psi ϭ 2.07 ϫ ϫ 10 6 Nm Ϫ2 is considered here.
Numerical calculations are carried out for a time-step Δτ ϭ ϭ 0.357 и 10 Ϫ4 s. Time variation of the deflection at the center of a simply supported plates subjected to a suddenly applied load is given in Fig. 5.
The peaks of the deflection and bending moment amplitudes are shifted to shorter time instants for the orthotropic homogeneous and laminated plates due to larger Young's moduli. They are largest for the laminated orthotropic plate. Since the mass density is the same in all plates, the elastic wave velocity is largest for the laminated plate. The maximum reduction of the deflection is achieved for the laminate plate, where the flexural rigidity is the largest.

Conclusion
A meshless local Petrov-Galerkin method is applied to laminate plates under mechanical loadings. The present computational method has no restriction on the number of the laminae and their material properties. The laminate plate bending problem is described by the Reissner-Mindlin theory. The Reissner-Mindlin theory reduces the original three-dimensional (3-D) thick plate problem to a 2-D problem. The weak-form on small subdomains with a Heaviside step function as the test functions is applied to derive local integral equations. After performing the spatial MLS approximation, a system of ordinary differential equations for certain nodal unknowns is obtained. Then, the system of the ordinary differential equations is solved by the Houbolt finite-difference scheme as a time-stepping method.
The proposed method is a truly meshless method, which requires neither domain elements nor background cells in either the interpolation or the integration. It is demonstrated numerically that the quality of the results obtained by the proposed MLPG method is very good. The degree of the agreement of our numerical results with those obtained by the FEM-ANSYS computer code ranges from good to excellent. Since in our illustrative examples only simple problems are analysed, only a regular node distribution has been used. However, a random location of nodes should be considered for further progress of the method.