Natural Frequencies of Shear Deformable Plates by Polyharmonic Splines

Abstract

We use the third-order shear deformation theory and a collocation technique with polyharmonic splines to predict natural frequencies of moderately thick isotropic plates. The natural frequencies of vibration are computed for various plates and compared with some available published results. Through numerical experiments, the capability and efficiency of the present method for eigenvalue problems are demonstrated, and the numerical accuracy and convergence are thoughtfully examined.

1. Introduction

Radial basis functions (RBFs) have recently proved to be an excellent technique for interpolating data and functions. A radial basis function 𝜙(‖𝑥−𝑥_𝑗‖) can be considered a spline that depends on the Euclidean distance between distinct data centers 𝑥_𝑗,𝑗=1,2,…,𝑁∈ℝ^𝑛, also called nodal or collocation points. Not only RBFs are adequate to scattered data approximation and in general to interpolation theory, as proposed by Hardy [1], but also excellent for solving partial differential equations (PDEs) as first introduced by Kansa [2].

Kansa proposed an unsymmetric RBF collocation method based upon multiquadric interpolation functions, in which the shape parameter is considered to be variable across the problem domain. The distribution of the shape parameter is obtained by an optimization approach, in which the value of the shape parameter is assumed to be proportional to the curvature of the unknown solution of the original partial differential equation. In this way, it is possible to reduce the condition number of the matrix at the expense of implementing an additional iterative algorithm. In the present work, we will implement the unsymmetric global collocation method in a form that is independent of this shape parameter, based on unshifted polyharmonic splines. In some respect this can be seen as a more stable form than multiquadrics.

Structures composed of laminated materials are among the most important structures used in modern engineering and, especially, in the aerospace industry. Such lightweight structures are also being increasingly used in civil, mechanical, and transportation engineering applications. The rapid increase of the industrial use of these structures has necessitated the development of new analytical and numerical tools that are suitable for the analysis and study of the mechanical behavior of such structures. The behavior of structures composed of advanced composite materials is considerably more complicated than for isotropic ones. The strong influences of anisotropy, the transverse stresses through the thickness of a laminate, and the stress distributions at interfaces are among the most important factors that affect the general performance of such structures. The use of shear deformation theories has been the topic of intensive research, as in [7–19], among many others.

The analysis of laminated plates by finite element methods is now considerably established. The use of alternative methods such as the meshless methods based on radial basis functions is attractive due to the absence of a mesh and the ease of collocation methods. The use of radial basis function for the analysis of structures and materials has been previously studied by numerous authors [20–31]. More recently the authors have applied RBFs to the static deformations of composite beams and plates [32–35]. Recent meshless techniques were discussed in [36–42].

In this paper the use of radial basis functions to isotropic and composite plates using a third-order shear deformation theory is investigated. The quality of the present method in predicting free vibrations of isotropic and laminated composite plates is compared and discussed with other methods in some numerical examples.

2. The Radial Basis Function Method

2.1. The Eigenproblem

Radial basis functions (RBF) approximations are grid-free numerical schemes that can exploit accurate representations of the boundary, are easy to implement, and can be spectrally accurate [43, 44].

In this section the formulation of a global unsymmetrical collocation RBF-based method to compute eigenvalues of elliptic operators is presented.

Consider a linear elliptic partial differential operator 𝐿 and a bounded region Ω in ℝ^𝑛 with some boundary 𝜕Ω. The eigenproblem looks for eigenvalues (𝜆) and eigenvectors (𝐮) that satisfy 𝐿𝐿𝐮+𝜆𝐮=0inΩ,_𝐵𝐮=0on𝜕Ω, (1) where 𝐿_𝐵 is a linear boundary operator. The eigenproblem of (1) is replaced by a finite-dimensional eigenvalue problem, based on RBF approximations.

The operator 𝐿 is approximated by a matrix that incorporates the boundary conditions and then solve the eigenvalues and eigenvectors of this matrix by standard techniques.

2.2. Radial Basis Functions

The radial basis function (𝜙) approximation of a function (𝐮) is given by 𝐮(𝐱)=^𝑁_𝑖=1𝛼_𝑖𝜙𝑥−𝑦_𝑖2,𝐱∈ℝ^𝑛, (2) where 𝑦_𝑖,𝑖=1,…,𝑁 is a finite set of distinct points (centers) in ℝ^𝑛. The coefficients 𝛼_𝑖 are calculated so that 𝐮 satisfies some boundary conditions. Although many other functions could be used, as illustrated by several authors [1, 2, 45], here we use unshifted polyharmonic splines in the form 𝜙(𝑟)=𝑟^2𝑚+1,𝑚∈𝑍, (3) where the Euclidean distance 𝑟 is real and nonnegative.

Considering 𝑁 distinct interpolations, and knowing 𝑢(𝑥_𝑗),𝑗=1,2,…,𝑁, we find 𝛼_𝑖 by the solution of a 𝑁×𝑁 linear system 𝐀𝛼=𝐮, where 𝐀=[𝜙(‖𝑥−𝑦_𝑖‖₂)]_𝑁×𝑁, 𝛼=[𝛼₁,𝛼₂,…,𝛼_𝑁]^𝑇 and 𝐮 = [𝑢(𝑥₁),𝑢(𝑥₂),…,𝑢(𝑥_𝑁)]^𝑇. The RBF interpolation matrix 𝐴 is positive definite for some RBFs [46] but in general provides ill-conditioned systems. This ill-conditioning problem is due to full matrices produced by the collocation method and worsens with the increase of number of nodes. However it was shown by Schaback that the more ill-conditioned the problem the better the solution, until a flatness of the function occurs. In our numerical experiments this ill-conditioning did not produce any inconvenience to quality of solution.

2.3. Solution of the Eigenproblem

We follow a simple scheme for the solution of the eigenproblem (1). We consider 𝑁_𝐼 nodes in the interior of the domain and 𝑁_𝐵 nodes on the boundary, with 𝑁=𝑁_𝐼+𝑁_𝐵.

We denote interpolation points by 𝑥_𝑖∈Ω,𝑖=1,…,𝑁_𝐼 and 𝑥_𝑖∈𝜕Ω,𝑖=𝑁_𝐼+1,…,𝑁. For the interior points we have that ^𝑁_𝑖=1𝛼_𝑖𝐿𝜙𝑥−𝑦_𝑖2=𝜆𝐮𝑥_𝑗,𝑗=1,2,…,𝑁_𝐼 (4) or 𝐿^𝐼𝛼=𝜆𝐮^𝐼, (5) where 𝐿^𝐼=𝐿𝜙𝑥−𝑦_{𝑖2𝑁𝐼×𝑁}, (6)

For the boundary conditions we have ^𝑁_𝑖=1𝛼_𝑖𝐿_𝐵𝜙𝑥−𝑦_𝑖2=0,𝑗=𝑁_𝐼+1,…,𝑁 (7) or 𝐁𝛼=0. (8)

Therefore we can write a finite-dimensional problem as a generalized eigenvalue problem 𝐿^𝐼𝐁𝛼𝐀=𝜆^𝐼𝟎𝛼, (9) where 𝐀^𝐼=𝜙𝑥_𝑁𝐼−𝑦_{𝑗2𝑁𝐼×𝑁},𝐁^𝐼=𝐿_𝐵𝜙𝑥_𝑁𝐼+1−𝑦_{𝑗2𝑁𝐵×𝑁}. (10)

3. Third-Order Theory

The displacement field for the third-order shear deformation theory of Hardy [1] is obtained as 𝑢𝑥,𝑦,𝑧,𝑡=𝑢₀𝑥,𝑦,𝑡+𝑧𝜃_𝑥−4𝑥,𝑦,𝑡3ℎ²𝑧³𝜃_𝑥+𝜕𝑤,𝑣𝜕𝑥𝑥,𝑦,𝑧,𝑡=𝑣₀𝑥,𝑦,𝑡+𝑧𝜃_𝑦−4𝑥,𝑦,𝑡3ℎ²𝑧³𝜃_𝑦+𝜕𝑤,𝑤𝜕𝑦𝑥,𝑦,𝑧,𝑡=𝑤₀𝑥,𝑦,𝑡, (11) where 𝑢 and 𝑣 are the in-plane displacements at any point (𝑥,𝑦,𝑧), 𝑢₀,𝑣₀ denote the in-plane displacement of the point (𝑥,𝑦,0) on the midplane, 𝑤 is the deflection, and 𝜃_𝑥 and 𝜃_𝑦 are the rotations of the normals to the midplane about the 𝑦- and 𝑥-axes, respectively.

The strain-displacement relationships are given as 𝜖_𝑥𝑥𝜖_𝑦𝑦𝛾_𝑥𝑦𝛾_𝑥𝑧𝛾_𝑦𝑧=𝜕𝑢𝜕𝑥𝜕𝑣𝜕𝑦𝜕𝑢+𝜕𝑦𝜕𝑣𝜕𝑥𝜕𝑢+𝜕𝑧𝜕𝑤𝜕𝑥𝜕𝑣+𝜕𝑧𝜕𝑤𝜕𝑦. (12)

Therefore strains can be expressed as 𝜖_𝑥𝑥𝜖_𝑦𝑦𝛾_𝑥𝑦=𝜖⁽⁰⁾_𝑥𝑥𝜖⁽⁰⁾_𝑦𝑦𝛾⁽⁰⁾_𝑥𝑦𝜖+𝑧𝜖^{(1)𝑥𝑥(1)}_𝑦𝑦𝛾⁽¹⁾_𝑥𝑦+𝑧³𝜖⁽³⁾_𝑥𝑥𝜖⁽³⁾_𝑦𝑦𝛾⁽³⁾_𝑥𝑦,𝛾_𝑥𝑧𝛾_𝑦𝑧=𝛾⁽⁰⁾_𝑥𝑧𝛾⁽⁰⁾_𝑦𝑧+𝑧²𝛾⁽²⁾_𝑥𝑧𝛾⁽²⁾_𝑦𝑧, (13) where 𝜖⁽⁰⁾_𝑥𝑥𝜖⁽⁰⁾_𝑦𝑦𝛾⁽⁰⁾_𝑥𝑦=𝜕𝑢₀𝜕𝑥𝜕𝑣₀𝜕𝑦𝜕𝑢₀+𝜕𝑦𝜕𝑣₀,𝜖𝜕𝑥_(1)𝑥𝑥𝜖⁽¹⁾_𝑦𝑦𝛾⁽¹⁾_𝑥𝑦=𝜕𝜃_𝑥𝜕𝑥𝜕𝜃_𝑦𝜕𝑦𝜕𝜃_𝑥+𝜕𝑦𝜕𝜃_𝑦,𝜖𝜕𝑥_(3)𝑥𝑥𝜖⁽³⁾_𝑦𝑦𝛾⁽³⁾_𝑥𝑦=−𝑐₁𝜕𝜃_𝑥+𝜕𝜕𝑥²𝑤₀𝜕𝑥²𝜕𝜃_𝑦+𝜕𝜕𝑦²𝑤₀𝜕𝑦²𝜕𝜃_𝑥+𝜕𝑦𝜕𝜃_𝑦𝜕𝜕𝑥+2²𝑤₀,𝛾𝜕𝑥𝜕𝑦_(0)𝑥𝑧𝛾⁽⁰⁾_𝑦𝑧=𝜕𝑤₀𝜕𝑥+𝜃_𝑥𝜕𝑤₀𝜕𝑦+𝜃_𝑦;𝛾⁽²⁾_𝑥𝑧𝛾⁽²⁾_𝑦𝑧=−𝑐₂𝜕𝑤₀𝜕𝑥+𝜃_𝑥𝜕𝑤₀𝜕𝑦+𝜃_𝑦 (14) and 𝑐₁=4/3ℎ²,𝑐₂=3𝑐₁.

A laminate can be manufactured from orthotropic layers (or plies) of preimpregnated unidirectional fibrous composite materials. Neglecting 𝜎_𝑧 for each layer, the stress-strain relations in the fiber local coordinate system can be expressed as 𝜎₁𝜎₂𝜏₁₂𝜏₂₃𝜏₃₁=𝑄₁₁𝑄₁₂𝑄000₁₂𝑄₂₂00000𝑄₃₃00000𝑄₄₄00000𝑄₅₅𝜀₁𝜀₂𝛾₁₂𝛾₂₃𝛾₃₁, (15) where subscripts 1 and 2 are, respectively, the fiber and the normal to fiber in-plane directions, 3 is the direction normal to the plate, and the reduced stiffness components, 𝑄_𝑖𝑗, are given by 𝑄₁₁=𝐸₁1−𝜈₁₂𝜈₂₁𝑄₂₂=𝐸₂1−𝜈₁₂𝜈₂₁𝑄₁₂=𝜈₂₁𝑄₁₁𝑄₃₃=𝐺₁₂𝑄₄₄=𝐺₂₃𝑄₅₅=𝐺₃₁𝜈₂₁=𝜈₁₂𝐸₂𝐸₁, (16) in which 𝐸₁,𝐸₂,𝜈₁₂,𝐺₁₂,𝐺₂₃, and 𝐺₃₁ are materials properties of the lamina.

By performing adequate coordinate transformation, the stress-strain relations in the global 𝑥-𝑦-𝑧 coordinate system can be obtained as 𝜎_𝑥𝑥𝜎_𝑦𝑦𝜏_𝑥𝑦𝜏_𝑦𝑧𝜏_𝑧𝑥=𝑄₁₁𝑄₁₂𝑄₁₆00𝑄₁₂𝑄₂₂𝑄₂₆00𝑄₁₆𝑄₂₆𝑄₆₆00000𝑄₄₄𝑄₄₅000𝑄₄₅𝑄₅₅𝜀_𝑥𝑥𝜀_𝑦𝑦𝛾_𝑥𝑦𝛾_𝑦𝑧𝛾_𝑧𝑥. (17)

The third-order theory of Hardy [15, 16] satisfies zero transverse shear stresses on the bounding planes.

The equations of motion of the third-order theory are derived from the principle of virtual displacements: 𝜕𝑁_𝑥𝑥+𝜕𝑥𝜕𝑁_𝑥𝑦𝜕𝑦=𝐼₀𝜕²𝑢₀𝜕𝑡²+𝐽₁𝜕²𝜃_𝑥𝜕𝑡²−𝑐₁𝐼₃𝜕²𝜕𝑡²𝜕𝑤₀,𝜕𝑥𝜕𝑁_𝑥𝑦+𝜕𝑥𝜕𝑁_𝑦𝑦𝜕𝑦=𝐼₀𝜕²𝑣₀𝜕𝑡²+𝐽₁𝜕²𝜃_𝑦𝜕𝑡²−𝑐₁𝐼₃𝜕²𝜕𝑡²𝜕𝑤₀,𝜕𝜕𝑦𝑄_𝑥+𝜕𝜕𝑥𝑄_𝑦𝜕𝑦+𝑐₁𝜕²𝑃_𝑥𝑥𝜕𝑥²𝜕+2²𝑃_𝑥𝑦+𝜕𝜕𝑥𝜕𝑦²𝑃_𝑦𝑦𝜕𝑦²+𝑞=𝐼₀𝜕²𝑤₀𝜕𝑡²−𝑐₁𝐼₆𝜕²𝜕𝑡²𝜕²𝑤₀𝜕𝑥²+𝜕²𝑤₀𝜕𝑦²+𝑐₁𝐼₃𝜕²𝜕𝑡²𝜕𝑢₀+𝜕𝑥𝜕𝑣₀𝜕𝑦+𝐽₄𝜕²𝜕𝑡²𝜕𝜃_𝑥+𝜕𝑥𝜕𝜃_𝑦,𝜕𝜕𝑦𝑀_𝑥𝑥+𝜕𝜕𝑥𝑀_𝑥𝑦−𝜕𝑦𝑄_𝑥=𝜕²𝜕𝑡²𝐽₁𝑢₀+𝐾₂𝜃_𝑥−𝑐₁𝐽₄𝜕𝑤₀,𝜕𝜕𝑥𝑀_𝑥𝑦+𝜕𝜕𝑥𝑀_𝑦𝑦−𝜕𝑦𝑄_𝑦=𝜕²𝜕𝑡²𝐽₁𝑣₀+𝐾₂𝜃_𝑦−𝑐₁𝐽₄𝜕𝑤₀,𝜕𝑦 (3) where 𝑞 is the external distributed load and with 𝑀_𝛼𝛽=𝑀_𝛼𝛽−𝑐₁𝑃_𝛼𝛽;𝑄_𝛼=𝑄_𝛼−𝑐₂𝑅_𝛼,𝐼_𝑖=^𝑁_𝑘=1𝑘^𝑘+1𝜌^(𝑘)(𝑧)^𝑖(),𝐽𝑑𝑧,𝑖=0,1,2,…,6_𝑖=𝐼_𝑖−𝑐₁𝐼_𝑖+2(),𝐾𝑖=1,4₂=𝐼₂−2𝑐₁𝐼₄+𝑐²₁𝐼₆,𝑐₁=43ℎ²,𝑐₂=4ℎ²=3𝑐₁, (19) where 𝛼,𝛽 take the symbols 𝑥,𝑦. The resultants (𝑁_𝑥𝑥,𝑁_𝑦𝑦,𝑁_𝑥𝑦) denote the in-plane force resultants, (𝑀_𝑥𝑥,𝑀_𝑦𝑦,𝑀_𝑥𝑦) the moment resultants, (𝑄_𝑥,𝑄_𝑦) the shear resultants and (𝑃_𝑥𝑥,𝑃_𝑦𝑦,𝑃_𝑥𝑦) and (𝑅_𝑥,𝑅_𝑦) denote the higher-order stress resultants: 𝑁_𝛼𝛽𝑀_𝛼𝛽𝑃_𝛼𝛽=^ℎ/2_−ℎ/2𝜎_𝛼𝛽1𝑧𝑧³𝑄𝑑𝑧,_𝛼𝑅_𝛼=^ℎ/2_−ℎ/2𝜎_𝛼𝑧1𝑧²𝑑𝑧, (20) where 𝛼,𝛽 take the symbols 𝑥,𝑦.

If needed, the first-order shear deformation theory equations are readily obtained from the third-order equations, just by putting 𝑐₁=0.

4. Numerical Examples

The eigenvalues presented in this problem are expressed in terms of the nondimensional frequency parameter 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2. Unlike first-order shear deformation theories the present approach does not incorporate shear correction factors in the constitutive matrix. The Poisson ratio, 𝜈, is taken to be 0.3 throughout the present problem. Figure 1 illustrates the grid scheme for 11×11 points.

Figure 1 

11×11 regular grid for discretization of the square laminated plate.

In Tables 1, 2, 3, 4, 5, and 6 we compare the present results with those of Liew et al. [3]. We consider here simply supported (SSSS) and clamped (CCCC) boundary conditions in all edges, thin and thick plates and square and rectangular plates. The results obtained in Tables 1 to 6 show that the present method presents very close results to those proposed by Liew et al. [3]. However, for thin plates we need more points than in thick plates in order to get convergent and accurate frequencies. The reason for this behaviour lies in the worse conditioning of thin plates. In Liew et al. [3] thin plates also showed slower convergence.

Table 1 

Convergence study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for rectangular isotropic plates with SSSS boundaries, ℎ/𝑎=0.001,𝑏/𝑎=1.0.

Table 2 

Convergence study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for rectangular isotropic plates with SSSS boundaries, ℎ/𝑎=0.2,𝑏/𝑎=1.0.

Table 3 

Convergence study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for rectangular isotropic plates with SSSS boundaries, ℎ/𝑎=0.001,𝑏/𝑎=2.5.

Table 4 

Convergence study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for rectangular isotropic plates with SSSS boundaries, ℎ/𝑎=0.2,𝑏/𝑎=2.5.

Table 5 

Convergence study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for rectangular isotropic plates with CCCC boundaries, ℎ/𝑎=0.001,𝑏/𝑎=1.0.

Table 6 

Convergence study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for rectangular isotropic plates with CCCC boundaries, ℎ/𝑎=0.2,𝑏/𝑎=1.0.

In Tables 7 and 8 we analyse thin and thick plates for SSSS and CCCC boundaries. Present results are compared with classical theory of Leissa [4], first-order theory of Dawe [5], and orthogonal polynomials of Liew et al. [3]. The present method shows very accurate results when compared to Liew et al. [3] and FSDT of Dawe [5], for both SSSS and CCCC boundary conditions.

Table 7 

Comparison study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for square isotropic plates with SSSS boundaries.

Table 8 

Comparison study of frequency parameters, 𝜆=(𝜔𝑏²/𝜋²)(𝜌𝑡/𝐷)^1/2 for square isotropic plates with CCCC boundaries, ℎ/𝑎=0.1.

5. Conclusions

Radial basis functions are increasingly popular in the analysis of science and engineering problems.

In this paper, we use the third-order shear deformation theory and a collocation technique with polyharmonic splines to predict natural frequencies of moderately thick isotropic plates. The natural frequencies of vibration are computed for various plates and compared with some available published results. The formulation for higher-order plates and their interpolation with polyharmonic splines has been presented.

Through numerical experiments, the capability and efficiency of the present method for eigenvalue problems are demonstrated, and the numerical accuracy and convergence are thoughtfully examined.

When compared to other RBFs, polyharmonic splines show good stability and accuracy. However, proper modelling of plates should consider both a good numerical technique and an adequate shear deformation theory.