A new meshless approach for bending analysis of thin plates with arbitrary shapes and boundary conditions

Abstract

An efficient and meshfree approach is proposed for the bending analysis of thin plates. The approach is based on the choice of a set of interior points, for each of which a basis function can be defined. Plate deflection is then approximated as the linear combination of those basis functions. Unlike traditional meshless methods, present basis functions are defined in the whole domain and satisfy the governing differential equation for plate. Therefore, no domain integration is needed, while the unknown coefficients of deflection expression could be determined through boundary conditions by using a collocation point method. Both efficiency and accuracy of the approach are shown through numerical results of plates with arbitrary shapes and boundary conditions under various loads.

Graphical abstract

Keywords

plate / bending / meshless method / collocation

1 Introduction

As a structural member, plate components are widely used in civil, mechanical, and aeronautical engineering due to light weight and high load-carrying capacity. To analyze the bending problems of plate, various methods have been developed, including both analytical and numerical approaches. Analytical methods are applicable only for plates with regular geometries and simple loads [1]. Therefore, analytical solutions are not available in most cases, and numerical approaches are required.

In the past decades, the finite difference method (FDM) [1], finite element method (FEM) [2] and boundary element method (BEM) [3] have been the most important numerical approaches for plate problems. Among them, FEM is commonly used in both research and engineering communities. FEM has attracted many scholars’ attention, and improvements have been made (for example, see Refs. [4–10]). When FEM is applied, discretization over the entire domain of a plate is necessary to generate a quality mesh, which is usually a time consuming process. To overcome the disadvantage of FEM, a meshless or meshfree method [11–15] has been developed recently. This method does not need any discretization of domain, and a local basis function is defined for appropriately chosen interior points. The approximated plate deflection may be constructed by using a technique such as moving least square. However, this process of construction is complex and time-consuming. As an improvement, an alternative form of meshfree method, called the line element-less method, has been introduced for the bending analysis of plates without any holes [16,17]. Several neoteric numerical methods, such as isogeometric analysis [18], numerical manifold method [19,20] and deep learning methods [21–23], have also been developed.

Presently, an efficient and meshless approach is proposed for the analysis of thin plates with arbitrary geometries and boundary conditions under various loads. The principle is based on selection of a set of points in the plate domain, and then for each point constructing a basis function. The basis function covers the whole domain and satisfies the governing differential equation of plate bending. Plate deflection is then approximated as the linear combination of those basis functions, while the unknown coefficients in the approximation can be determined directly through boundary conditions by using a collocation point method. Since the approximated deflection satisfies the governing differential equation exactly, no domain integration is needed, which will largely reduce calculation.

The remaining of this paper is organized as follows: Section 2 presents the governing differential equation of plate bending, as well as the expression of internal forces. Section 3 derives the formulation of the present method in detail. In Section 4, boundary collocation point method is introduced to determine unknown coefficients. Section 5 shows numerical results and discussions. Section 6 makes a conclusion.

2 Governing differential equations of plate bending

Consider the bending of a homogeneous isotropic thin plate. The governing differential equation for deflection of plate in Cartesian coordinates is

(1)D{\mathrm{\nabla }}^2{\mathrm{\nabla }}^{2}w=D\left(\frac{{\mathrm{\partial }}^{4}w}{\mathrm{\partial }{x}^4}+2\frac{{\mathrm{\partial }}^{4}w}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^{4}w}{\mathrm{\partial }{y}^4}\right)=q(x,y),

where D=E{h}^3/12(1-{\mu }^2) and is the flexural rigidity of the plate, determined by the elastic modulus E, the plate thickness h, and the Poisson’s ratio \mu. {\mathrm{\nabla }}^2={\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}}+{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}} and is the Laplacian operator. q(x,y) is the transverse load. The internal forces can be written in terms of the deflection w as follows:

(2)\begin{array}{rl}& M_x=-D\left(\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}+\mu \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}\right),\\ & M_y=-D\left(\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}+\mu \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}\right),\\ & M_{xy}=M_{yx}=-D(1-\mu )\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }x\mathrm{\partial }y},\\ & Q_x=-D\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\mathrm{\nabla }}^{2}w\right),\\ & Q_y=-D\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\mathrm{\nabla }}^{2}w\right).\end{array}

To facilitate derivation, it is convenient to express the governing differential equation and internal forces in polar coordinates. The Laplacian operator takes the form

(3){\mathrm{\nabla }}^{2}w=\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}=\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{r}^2}+\frac{1}{r}\frac{\mathrm{\partial }w}{\mathrm{\partial }r}+\frac{1}{r^2}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{\theta }^2}.

The governing differential equation Eq. (1) becomes

(4)\begin{array}{rl}D{\mathrm{\nabla }}^2{\mathrm{\nabla }}^{2}w=& D\left(\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{r}^2}+\frac{1}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}+\frac{1}{r^2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\theta }^2}\right)\\ & \left(\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{r}^2}+\frac{1}{r}\frac{\mathrm{\partial }w}{\mathrm{\partial }r}+\frac{1}{r^2}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{\theta }^2}\right)=q(r,\theta ).\end{array}

The moments, twisting moments and shear forces can be written as

(5)\begin{array}{rl}& M_{\mathrm{r}}=-D\left[\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{r}^2}+\mu \left(\frac{1}{r}\frac{\mathrm{\partial }w}{\mathrm{\partial }r}+\frac{1}{r^2}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{\theta }^2}\right)\right],\\ & M_{\text{θ}}=-D\left[\left(\frac{1}{r}\frac{\mathrm{\partial }w}{\mathrm{\partial }r}+\frac{1}{r^2}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{\theta }^2}\right)+\mu \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{r}^2}\right],\\ & M_{\mathrm{r}\text{θ}}=M_{\text{θ}\mathrm{r}}=-D(1-\mu )\left[\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{1}{r}\frac{\mathrm{\partial }w}{\mathrm{\partial }\theta }\right)\right],\\ & Q_{\mathrm{r}}=-D\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left({\mathrm{\nabla }}^{2}w\right),\\ & Q_{\text{θ}}=-D\frac{1}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\left({\mathrm{\nabla }}^{2}w\right).\end{array}

3 Formulation of the present method

According to the theory of differential equations, the general solution for the governing differential equation Eq. (1) can be obtained by the superposition of a particular solution of the equation and the general solution of a corresponding homogeneous equation.

Firstly, a particular solution is obtained due to a concentrated force. Referring to Fig. 1, assume that a concentrated force P acts at an arbitrary point A_0(x_0,y_0). In this case, the differential equation Eq. (1) can be written as

(6)D{\mathrm{\nabla }}^2{\mathrm{\nabla }}^{2}w=P\cdot \delta (A,A_0),

where \delta (A,A_0) is the two-dimensional \delta-function, and

(7)\delta (A,A_0)=0(A\ne A_0),\iint \delta (A,A_0)\mathrm{d}S=1.

where dS = dxdy. The particular solution of Eq. (6) could take the following form [12]

(8){\tilde{w}}_1\left(x,y\right)=\frac{P}{8\text{π}D}{R}^2\ln R,

where R^2=(x-x_0{)}^2+(y-y_0{)}^2.

Using the above equation, the particular solution due to generally distributed loading q(x_0,y_0) is then obtained, as

(9)\begin{array}{rl}{\tilde{w}}_2\left(x,y\right)=& \underset{S}{\iint }\{\frac{q(x_0,y_0)}{8\text{π}D}\left[{(x-x_0)}^2+{(y-y_0)}^2\right]\\ & \ln \left(\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}\right)\}\mathrm{d}{x}_0\mathrm{d}{y}_0.\end{array}

Specially, for uniform loads, the particular solution could take the simple form

(10){\tilde{w}}_3\left(x,y\right)=\frac{q_0}{64D}{\left(x^2+y^2\right)}^2.

Next, the homogeneous equation D{\mathrm{\nabla }}^2{\mathrm{\nabla }}^{2}w=0 is to be discussed. Its solution can take the form [1]

(11)\begin{array}{rl}{w}_0\left(r,\theta \right)=& \left(a_0+a_0^{\ast }\ln r+b_0{r}^2+b_0^{\ast }{r}^2\ln r\right)\\ & +\left(a_{1}r+a_1^{\ast }{r}^{-1}+b_1{r}^3+b_1^{\ast }r\ln r\right)\cos \theta \\ & +\sum _{k=2}^n\left(a_k{r}^k+a_k^{\ast }{r}^{-k}+b_k{r}^{2+k}+b_k^{\ast }{r}^{2-k}\right)\cos k\theta \\ & +\left(c_{1}r+c_1^{\ast }{r}^{-1}+d_1{r}^3+d_1^{\ast }r\ln r\right)\sin \theta \\ & +\sum _{k=2}^n\left(c_k{r}^k+c_k^{\ast }{r}^{-k}+d_k{r}^{2+k}+d_k^{\ast }{r}^{2-k}\right)\sin k\theta ,\end{array}

where a_k,a_k^{\ast },b_k,b_k^{\ast },c_k,c_k^{\ast },d_k,d_k^{\ast } are unknown coefficients.

If the origin of the coordinate system is located in the domain of the plate, Eq. (11) can be simplified further. In this case, the deflection w_0, the slope \mathrm{\partial }{w}_0/\mathrm{\partial }r and the moment should be finite at r=0. Moreover, the deflection w_0 is independent of \theta when r=0. Thus, a_0^{\ast }=a_k^{\ast }=b_k^{\ast }=c_k^{\ast }=d_k^{\ast }=0(k⩾1), and Eq. (11) becomes

(12)\begin{array}{rl}{w}_0\left(r,\theta \right)=& a_0+b_0{r}^2+b_0^{\ast }{r}^2\ln r\\ & +\sum _{k=1}^n\left(a_k{r}^k+b_k{r}^{2+k}\right)\cos k\theta \\ & +\sum _{k=1}^n\left(c_k{r}^k+d_k{r}^{2+k}\right)\sin k\theta .\end{array}

Note that w_0(r,\theta ) only represents the effect of boundary on deflection and should contain no effect of loading. Therefore, when a closed-circuit integral is performed for the lateral shear force in Eq. (5) along a small circle around r=0, there should be

(13){\displaystyle \oint Q_{\text{r}}\text{d}s=-{\displaystyle {\int }_0^{2\text{π}}D\frac{\mathrm{\partial }}{\mathrm{\partial }r}}}\left({\mathrm{\nabla }}^2{w}_0\right)\cdot r\text{d}\theta =-8\text{π}D{b}_0^{\ast }=0.

Thus, b_0^{\ast }=0. The solution (12) is simplified as

(14)w_0\left(r,\theta \right)=\sum _{k=0}^n\left(a_k{r}^k+b_k{r}^{2+k}\right)\cos k\theta +\sum _{k=1}^n\left(c_k{r}^k+d_k{r}^{2+k}\right)\sin k\theta .

The above deflection representation in polar coordinates could be converted to Cartesian form. According to the theory of complex analysis, harmonic polynomials are defined as follows

(15)\begin{array}{rl}& P_k(x,y)=r^k\mathrm{c}\mathrm{o}\mathrm{s}k\theta =\mathrm{R}\mathrm{e}{(x+\text{i}y)}^k;\\ & Q_k(x,y)=r^k\mathrm{s}\mathrm{i}\mathrm{n}k\theta =\mathrm{I}\mathrm{m}{(x+\text{i}y)}^k,\end{array}

where ‘Re’ and ‘Im’ represent the real and imaginary part of a complex variable, respectively. They have the following recursive form

(16)\begin{array}{rl}& P_k(x,y)=x{P}_{k-1}-y{Q}_{k-1};\\ & Q_k(x,y)=y{P}_{k-1}+x{Q}_{k-1}(\mathrm{\forall }k>0)\end{array},

where P_0=1, Q_0=0, and P_1=x, Q_1=y. Their derivatives can also be expressed recursively as

(17)\begin{array}{rl}& \frac{\mathrm{\partial }{P}_k}{\mathrm{\partial }x}=k{P}_{k-1},\frac{\mathrm{\partial }{P}_k}{\mathrm{\partial }y}=-k{Q}_{k-1},\\ & \frac{\mathrm{\partial }{Q}_k}{\mathrm{\partial }x}=k{Q}_{k-1},\frac{\mathrm{\partial }{Q}_k}{\mathrm{\partial }y}=k{P}_{k-1}\left(\mathrm{\forall }k>0\right)\end{array}.

Substituting Eq. (15) into Eq. (14), w_0(r,\theta ) could be recast as

(18)\begin{array}{rl}{w}_0(x,y)=& \sum _{k=0}^n\left[a_k{P}_k(x,y)+b_k\left(x^2+y^2\right)P_k(x,y)\right]\\ & +\sum _{k=1}^n\left[c_k{Q}_k(x,y)+d_k\left(x^2+y^2\right)Q_k(x,y)\right].\end{array}

Clearly, \mathrm{\forall }{x}_j,y_j( j=1,2,\dots ,m, as shown in Fig. 2), w_j(x,y)=w_0(x-x_j,y-y_j) satisfies the biharmonic equation D{\mathrm{\nabla }}^2{\mathrm{\nabla }}^{2}w=0. Hence, w_j(x,y)( j=1,2,\dots ,m) can be chosen as the basis functions. The general solution of D{\mathrm{\nabla }}^2{\mathrm{\nabla }}^{2}w=0 is expressed as

(19)\overline{w}(x,y)=\sum _{j=1}^m{w}_j(x,y)=\sum _{j=1}^m{w}_0(x-x_j,y-y_j),

where m denotes the number of interior points.

For convenience, assume that X_j=x-x_j, Y_j=y-y_j. w_j(x,y) can be written as

(20)\begin{array}{rl}{w}_j(x,y)=& a_{0j}+b_{0j}\left(X_j^2+Y_j^2\right)\\ & +\sum _{k=1}^n\left[a_{kj}{P}_k\left(X_j,Y_j\right)+b_{kj}\left(X_j^2+Y_j^2\right)P_k\left(X_j,Y_j\right)\right]\\ & +\sum _{k=1}^n\left[c_{kj}{Q}_k\left(X_j,Y_j\right)+d_{kj}\left(X_j^2+Y_j^2\right)Q_k\left(X_j,Y_j\right)\right]\end{array}.

Substituting the above equation into Eq. (19), we obtain

(21)\begin{array}{rl}\overline{w}(x,y)=& a_0+\sum _{j=1}^m{b}_{0j}\left(X_j^2+Y_j^2\right)\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[a_{kj}{P}_k\left(X_j,Y_j\right)+b_{kj}\left(X_j^2+Y_j^2\right)P_k\left(X_j,Y_j\right)\right]\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[c_{kj}{Q}_k\left(X_j,Y_j\right)+d_{kj}\left(X_j^2+Y_j^2\right)Q_k\left(X_j,Y_j\right)\right],\end{array}

where a_0=a_{01}+a_{02}+\cdots +a_{\text{0m}}.

Finally, the plate deflection can be represented as

(22)w(x,y)=\tilde{w}(x,y)+\overline{w}(x,y),

where \tilde{w}(x,y) is given in Eqs. (8), (9), or (10) for differently distributed loads. The coefficients a_0,b_{0j},a_{kj},b_{kj},c_{kj},d_{kj} in Eq. (22) are unknown and remain to be determined through boundary conditions.

The simplified form Eq. (14) of general solution Eq. (11) requires that the plate domain includes the origin r=0. However, the introduction of interior points lifts the constraint, as these interior points play the role of local origins. As a result, the final general deflection expression Eq. (22) is independent of the choice of coordinate system.

It should be noted that Eq. (22) is quite different from the traditional meshless method in that the present approximation of deflection satisfies the governing differential equation Eq. (1) exactly. The unknown coefficients in Eq. (22) could be determined through boundary conditions only, and therefore no domain integration is needed, which will reduce the cost of calculation hugely.

Further, introducing the following expressions

(23)\begin{array}{rl}& {\mathit{N}}_j(x,y)=\left\{\begin{array}{c}{X}_j^2+Y_j^2\\ ⋮\\ P_k(X_j,Y_j)\\ (X_j^2+Y_j^2)P_k(X_j,Y_j)\\ Q_k(X_j,Y_j)\\ (X_j^2+Y_j^2)Q_k(X_j,Y_j)\\ ⋮\end{array}\right\},\\ & {\mathit{\delta }}_{\mathbf{j}}=\left\{\begin{array}{c}{b}_{0j}\\ ⋮\\ a_{kj}\\ b_{kj}\\ c_{kj}\\ d_{kj}\\ ⋮\end{array}\right\}\left(\begin{array}{rl}& k=1,2,\cdots ,n\\ & j=1,2,\cdots ,m\end{array}\right).\end{array}

Equation (22) could be rewritten in compact matrix form as

(24)w\left(x,y\right)=\tilde{w}\left(x,y\right)+\left[\mathit{N}\right]\left\{\mathit{\delta }\right\},

where [\mathit{N}]=[1,{\mathit{N}}_1^{\mathrm{T}},{\mathit{N}}_2^{\mathrm{T}},\cdots ,{\mathit{N}}_m^{\mathrm{T}}], \{\mathit{\delta }\}=[a_0,{\mathit{\delta }}_1^{\mathrm{T}},{\mathit{\delta }}_2^{\mathrm{T}},\cdots ,{\mathit{\delta }}_m^{\mathrm{T}}{]}^{\mathrm{T}}.

4 Boundary collocation point

To determine the unknown coefficients in Eq. (22), boundary conditions should be considered. At present, the boundary conditions are imposed by the collocation point method. Assume that the residual along the plate boundary is R(s). Choose N discrete boundary points s_{\gamma }(\gamma =1,2,\dots ,N) and let R(s_{\gamma })=0, which leads to a system of linear algebraic equations for the unknown coefficients.

For a generally curvilinear edge, suppose that n and t are the outward unit normal and tangent vectors at a point on the edge, and \alpha is the angle between n and the x axis, resulting in the following boundary conditions.

(i) Clamped edge

(25)w=0,\frac{\mathrm{\partial }w}{\mathrm{\partial }n}=\cos \alpha \frac{\mathrm{\partial }w}{\mathrm{\partial }x}+\sin \alpha \frac{\mathrm{\partial }w}{\mathrm{\partial }y}=0.

These equations are imposed at N_1 discrete points s_{\gamma }(x_{\gamma },y_{\gamma }) on the edge, leading to a system of linear algebraic equations

(26)w\left(x_{\gamma },y_{\gamma }\right)=0,\frac{\mathrm{\partial }w(x_{\gamma },y_{\gamma })}{\mathrm{\partial }n}=0\left(\gamma =1,2,\dots ,N_1\right).

(ii) Simply-supported edge

(27)w=0,M_n=0,

where M_n(x,y) denotes the bending moment at normal direction

(28)M_n=M_x{\cos }^2\alpha +M_y{\sin }^2\alpha +M_{xy}\sin 2\alpha .

Substituting Eq. (2) into Eq. (28), M_n(x,y) could be expressed as

(29)M_n=-D\left(I_1\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}+I_2\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}+I_3\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }x\mathrm{\partial }y}\right),

where

(30)\begin{array}{rl}& I_1=\mu +\left(1-\mu \right){\cos }^2\alpha ,\\ & I_2=\mu +\left(1-\mu \right){\sin }^2\alpha ,\\ & I_3=\left(1-\mu \right)\sin 2\alpha .\end{array}

The boundary conditions at N_2 discrete points s_{\gamma }(x_{\gamma },y_{\gamma }) have the form

(31)w\left(x_{\gamma },y_{\gamma }\right)=0,M_n\left(x_{\gamma },y_{\gamma }\right)=0\left(\gamma =1,2,\dots ,N_2\right).

(iii) Free edge

(32)M_n=0,V_n=Q_n+\frac{\mathrm{\partial }{M}_{nt}}{\mathrm{\partial }s}=0,

where V_n(x,y) is the effective shear force, M_{nt}(x,y) and Q_n(x,y) represent the twisting moment and the shearing force

(33)\begin{array}{rl}& M_{nt}=\left(M_y-M_x\right)\cos \alpha \sin \alpha +M_{xy}\cos 2\alpha ,\\ & Q_n=Q_x\cos \alpha +Q_y\sin \alpha .\end{array}

Introducing Eq. (2)

(34)V_n=-D\left(H_1\frac{{\mathrm{\partial }}^{3}w}{\mathrm{\partial }{x}^3}+H_2\frac{{\mathrm{\partial }}^{3}w}{\mathrm{\partial }{y}^3}+H_3\frac{{\mathrm{\partial }}^{3}w}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+H_4\frac{{\mathrm{\partial }}^{3}w}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right),

where

(35)\begin{array}{rl}& H_1=\cos \alpha +\left(1-\mu \right)\cos \alpha {\sin }^2\alpha ,\\ & H_3=\cos \alpha +\left(1-\mu \right)({\cos }^3\alpha -2\cos \alpha {\sin }^2\alpha ),\\ & H_2=\sin \alpha +\left(1-\mu \right){\cos }^2\alpha \sin \alpha ,\\ & H_4=\sin \alpha +\left(1-\mu \right)({\sin }^3\alpha -2{\cos }^2\alpha \sin \alpha ).\end{array}\}

The boundary conditions at N_3 discrete points s_{\gamma }(x_{\gamma },y_{\gamma }) become

(36)M_n\left(x_{\gamma },y_{\gamma }\right)=0,V_n\left(x_{\gamma },y_{\gamma }\right)=0\left(\gamma =1,2,\dots ,N_3\right).

In the above equations, N_1+N_2+N_3=N. By solving the linear algebraic system of Eqs. (26), (31), and (36) together, the unknown coefficients are obtained. Further substitution in Eqs. (22) and (2) leads to the complete solution of plate deflection and internal forces.

5 Numerical examples

This section proposed the approach to analyze a few benchmark examples. The obtained results are compared with exact solutions or FEM numerical solutions given by ABAQUS. Remarkably, the proposed approach demonstrates excellent computational efficiency as its computational time is less than one-tenth that of FEM by ABAQUS with these benchmark tests.

5.1 Circular plate subjected to an eccentric concentrated load

As shown in Fig. 3, a circular plate with clamped edge is subjected to a unit concentrated load P=1 at point A_0(a/2,0). The radius of the plate is r=a.

The boundary conditions referring to Fig. 3 are represented by

(37){\left(w\right)}_{r=a}=0,{\left(\mathrm{\partial }w/\mathrm{\partial }n\right)}_{r=a}=0.

In this case x_0=a/2, y_0=0. According to Eqs. (8), (21), and (22), the deflection has the form

(38)\begin{array}{rl}w(x,y)=& \frac{P}{8\text{π}D}{R}^2\ln R+a_0+\sum _{j=1}^m{b}_{0j}\left(X_j^2+Y_j^2\right)\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[a_{kj}{P}_k\left(X_j,Y_j\right)+b_{kj}\left(X_j^2+Y_j^2\right)P_k\left(X_j,Y_j\right)\right]\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[c_{kj}{Q}_k\left(X_j,Y_j\right)+d_{kj}\left(X_j^2+Y_j^2\right)Q_k\left(X_j,Y_j\right)\right],\end{array}

where R^2=(x-a/2{)}^2+y^2.

We choose n=5, m=1, and (x_j,y_j)=(0,0). Following the procedure described in Section 4, the unknown coefficients can be obtained. When a=1, the exact solution [24] of this problem is

(39)\begin{array}{rl}{w}^{\ast }(r,\theta )=& \frac{1}{16\text{π}D}[\left(r^2+\frac{1}{4}-r\cos \theta \right)\\ & \ln \frac{r^2+1/4-r\cos \theta }{1+r^2/4-r\cos \theta }+\frac{3}{4}\left(1-r^2\right)].\end{array}

Figure 4 shows the contours of the deflected surface. This figure clearly shows the variation of deflection in the whole domain, with the maximum deflection appearing near the loading point. Figure 5 exhibits the deflection profile at y=0. In comparison with the analytical solution, good agreement of present results can be observed.

Figure 6 checks the convergence rate of the present method. Let {\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}} denote the relative error of the maximum deflection. It can be seen that with the increase of both terms n and the number of interior points m, present results converge rapidly to the exact solution. Specifically, for a definite m, the absolute value of {\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}} decreases almost monotonically with the increase of n, and becomes negligible when n is greater than 5. For a definite n greater than 3, the absolute value of {\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}} decreases monotonically with the increase of m, and becomes negligible when m is greater than 3.

Figures 7 and 8 show deflection variation along the boundary. Since the deflection on the boundary should be zero, the two figures represent the absolute errors of deflection on the boundary, reflecting both the convergence and accuracy of present results. Again, the absolute errors are negligible. The present method can achieve good convergence and accuracy, even with very few terms and interior points.

5.2 Sectorial plate subjected to uniformly distributed load

A sectorial plate (Fig. 9) with radius r=a has two clamped straight edges \theta =\pm \text{π}/4 and a free circular arc edge. Suppose that a uniformly distributed load q_0=1 is applied, and the Poisson’s ratio \mu =0.3.

The boundary conditions are

(40){\left(w\right)}_{\theta =\pm \pi /4}=0,{\left(\mathrm{\partial }w/\mathrm{\partial }n\right)}_{\theta =\pm \pi /4}=0,

(41){\left(M_n\right)}_{r=a}=0,{\left(Q_n+\frac{\mathrm{\partial }{M}_{nt}}{\mathrm{\partial }s}\right)}_{r=a}=0.

According to Eqs. (10), (21), and (22), the deflection is assumed to be

(42)\begin{array}{rl}w(x,y)=& \frac{q_0}{64D}{\left(x^2+y^2\right)}^2+a_0+\sum _{j=1}^m{b}_{0j}\left(X_j^2+Y_j^2\right)\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[a_{kj}{P}_k\left(X_j,Y_j\right)+b_{kj}\left(X_j^2+Y_j^2\right)P_k\left(X_j,Y_j\right)\right]\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[c_{kj}{Q}_k\left(X_j,Y_j\right)+d_{kj}\left(X_j^2+Y_j^2\right)Q_k\left(X_j,Y_j\right)\right].\end{array}

The unknown coefficients are obtained with n=20, m=1 and (x_j,y_j)=(0,0). Figure 10 shows the contours of deflected surface. In this case, the maximum deflection occurs at the middle point of the circular arc edge. From the deflection profile at y=0 (Fig. 11), excellent agreement is observed with the FEM solution obtained by ABAQUS (with 128 quadratic quadrilateral plate elements and 421 nodes). The relative error of maximum deflection is only -0.15\mathrm{\%}.

5.3 Cantilever rectangular plate subjected to concentrated loads

A cantilever rectangular plate is shown in Fig. 12. Suppose that the plate is subjected to a pair of unit concentrated loads P=1 at the points (\pm 3a/8,3b/4) symmetrically or anti-symmetrically. The Poisson’s ratio \mu =0.3.

The boundary conditions for the clamped edge and free edges are

(43){\left(w\right)}_{y=0}=0;{\left(\mathrm{\partial }w/\mathrm{\partial }y\right)}_{y=0}=0(-a/2⩽x⩽a/2),

(44){\left(M_x\right)}_{x=\pm a/2}=0,{\left(Q_x+\frac{\mathrm{\partial }{M}_{xy}}{\mathrm{\partial }y}\right)}_{x=\pm a/2}=0,(0⩽y⩽b)

(45){\left(M_y\right)}_{y=b}=0,{\left(Q_y+\frac{\mathrm{\partial }{M}_{yx}}{\mathrm{\partial }x}\right)}_{y=b}=0.(-a/2⩽x⩽a/2)

At two corners (\pm a/2,b), it needs to satisfy

(46){\left({\mathrm{\partial }}^{2}w/\mathrm{\partial }x\mathrm{\partial }y\right)}_{x=\pm a/2,y=b}=0.

The deflection has the following form

(47)\begin{array}{rl}w(x,y)=& \frac{P}{8\text{π}D}{R}_1^2\ln R_1\pm \frac{P}{8\text{π}D}{R}_2^2\ln R_2+a_0\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[a_{kj}{P}_k\left(X_j,Y_j\right)+b_{kj}\left(X_j^2+Y_j^2\right)P_k\left(X_j,Y_j\right)\right]\\ & +\sum _{j=1}^m{b}_{0j}\left(X_j^2+Y_j^2\right)\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[c_{kj}{Q}_k\left(X_j,Y_j\right)+d_{kj}\left(X_j^2+Y_j^2\right)Q_k\left(X_j,Y_j\right)\right],\end{array}

where R_1^2=(x-3a/8{)}^2+(y-3b/4{)}^2, R_2^2=(x+3a/8{)}^2+ (y-3b/4{)}^2. Consider the case of a=b. Choose four interior points located at (\pm 0.2a,0.3a) and (\pm 0.2a,0.7a). The unknown coefficients are obtained with n=16 and m=4.

Figures 13 and 14 show the deflected surfaces of symmetric and antisymmetric deformations, respectively. It is observed that the maximum deflections in both cases appear at the free corners. The deflection profile at x=0 for symmetric deformation (Fig. 15) and x=a/2 for antisymmetric deformation (Fig. 16) are in good agreement with the FEM solution obtained by ABAQUS (with 256 quadratic quadrilateral plate elements and 833 nodes). The absolute value of relative errors for maximum deflections are all within 0.5\mathrm{\%}.

5.4 Cantilever triangular plates subjected to concentrated load

Shown in Fig. 17 is a cantilever isosceles triangular plate with clamped bottom edge subjected to a unit lateral concentrated load P=1 at the point (0,b/2). \mu =0.3.

The boundary conditions are

(48){\left(w\right)}_{y=0}=0,{\left(\mathrm{\partial }w/\mathrm{\partial }y\right)}_{y=0}=0,(-a⩽x⩽a)

(49){\left(M_n\right)}_{x\mathrm{c}\mathrm{o}\mathrm{t}\varphi +y=b}=0,{\left(Q_n+\frac{\mathrm{\partial }{M}_{nt}}{\mathrm{\partial }s}\right)}_{x\cot \varphi +y=b}=0.(-a⩽x⩽a)

At the corner (0,b),

(50){\left[\left(I_1-I_2\right)\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{I_3}{2}\left(\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}-\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}\right)\right]}_{x=0,y=b}=0.

The deflection expression takes the same form as the circular plate, as shown in Eq. (38), where R^2=x^2+ (y-b/2{)}^2. Further, assume that a=b. Choose (0,a/3) and (0,2a/3) for (x_j,y_j). The unknown coefficients are obtained with n=12 and m=2.

Figure 18 shows the deflected surface. The maximum deflection occurs at the corner. Again, the deflection profile along the symmetry axis x=0 agrees well (Fig. 19) with the solution of FEM by ABAQUS (168 quadratic quadrilateral plate elements and 549 nodes). The relative error of maximum deflection is -0.29\mathrm{\%}.

5.5 Nonhomogeneous cantilever rectangular plate under uniform load

A nonhomogeneous cantilever rectangular plate (Fig. 20) is subjected to a uniform pressure of intensity q_0=1. The ratio of flexural rigidity between the two domains is D_2/D_1=0.512. The Poisson’s ratio \mu =0.3.

Referring to Fig. 20, the boundary conditions for domain I are as follows:

(51){\left(w_1\right)}_{y=0}=0,{\left(\mathrm{\partial }{w}_1/\mathrm{\partial }y\right)}_{y=0}=0(-a/2⩽x⩽a/2),

(52){\left(M_x^{\mathrm{I}}\right)}_{x=\pm a/2}=0;{\left(Q_x^{\mathrm{I}}+\frac{\mathrm{\partial }{M}_{xy}^{\mathrm{I}}}{\mathrm{\partial }y}\right)}_{x=\pm a/2}=0.(0⩽y⩽b/2)

The boundary conditions for domain II can be written as

(53){\left(M_x^{\mathrm{I}\mathrm{I}}\right)}_{x=\pm a/2}=0;{\left(Q_x^{\mathrm{I}\mathrm{I}}+\frac{\mathrm{\partial }{M}_{xy}^{\mathrm{I}\mathrm{I}}}{\mathrm{\partial }y}\right)}_{x=\pm a/2}=0,(b/2⩽y⩽b)

(54){\left(M_y^{\mathrm{I}\mathrm{I}}\right)}_{y=b}=0;{\left(Q_y^{\mathrm{I}\mathrm{I}}+\frac{\mathrm{\partial }{M}_{yx}^{\mathrm{I}\mathrm{I}}}{\mathrm{\partial }x}\right)}_{y=b}=0.(-a/2⩽x⩽a/2)

At two corners (\pm a/2,b),

(55){\left(\frac{{\mathrm{\partial }}^2{w}_2}{\mathrm{\partial }x\mathrm{\partial }y}\right)}_{x=\pm a/2,y=b}=0.

The connection conditions along the line y=b/2(-a/2⩽x⩽a/2) are

(56)\begin{array}{rl}& {\left(w_1-w_2\right)}_{y=\frac{b}{2}}=0,\\ & {\left(\frac{\mathrm{\partial }{w}_1}{\mathrm{\partial }y}-\frac{\mathrm{\partial }{w}_2}{\mathrm{\partial }y}\right)}_{y=b/2}=0\\ & {\left(M_y^{\mathrm{I}}-M_y^{\mathrm{I}\mathrm{I}}\right)}_{y=\frac{b}{2}}=0,\\ & {\left(M_{yx}^{\mathrm{I}}-M_{yx}^{\mathrm{I}\mathrm{I}}\right)}_{y=\frac{b}{2}}=0,\\ & {\left(Q_y^{\mathrm{I}}-Q_y^{\mathrm{I}\mathrm{I}}\right)}_{y=\frac{b}{2}}=0.\end{array}\}

In domain I, the deflection w_1 takes the form

(57)\begin{array}{rl}{w}_1(x,y)=& \frac{q_0}{64D_1}{\left(x^2+y^2\right)}^2+a_0\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[a_{kj}{P}_k\left(X_j,Y_j\right)+b_{kj}\left(X_j^2+Y_j^2\right)P_k\left(X_j,Y_j\right)\right]\\ & +\sum _{j=1}^m{b}_{0j}\left(X_j^2+Y_j^2\right)\\ & +\sum _{j=1}^m\sum _{k=1}^n\left[c_{kj}{Q}_k\left(X_j,Y_j\right)+d_{kj}\left(X_j^2+Y_j^2\right)Q_k\left(X_j,Y_j\right)\right],\end{array}

where X_j=x-x_j, Y_j=y-y_j.

In domain II, the deflection w_2 is

(58)\begin{array}{rl}{w}_2(x,y)=& \frac{q_0}{64D_2}{\left(x^2+y^2\right)}^2+a_0^{\ast }\\ & +\sum _{\beta =1}^m\sum _{k=1}^n\left[a_{k\beta }^{\ast }{P}_k\left(X_{\beta },Y_{\beta }\right)+b_{k\beta }^{\ast }\left(X_{\beta }^2+Y_{\beta }^2\right)P_k\left(X_{\beta },Y_{\beta }\right)\right]\\ & +\sum _{\beta =1}^m{b}_{0\beta }^{\ast }\left(X_{\beta }^2+Y_{\beta }^2\right)\\ & +\sum _{\beta =1}^m\sum _{k=1}^n\left[c_{k\beta }^{\ast }{Q}_k\left(X_{\beta },Y_{\beta }\right)+d_{k\beta }^{\ast }\left(X_{\beta }^2+Y_{\beta }^2\right)Q_k\left(X_{\beta },Y_{\beta }\right)\right],\end{array}

where X_{\beta }=x-x_{\beta }, Y_{\beta }=y-y_{\beta }.

Assume that a=b and D_1=D. Choose interior points (\pm a/4,a/4) and (0,a/4) for (x_j,y_j), (\pm a/4,3a/4) and (0,3a/4) for (x_{\beta },y_{\beta }). The unknown coefficients are obtained with n=20 and m=3.

Figure 21 shows the deflected surface, with the maximum deflection at the middle point of free edge y=b. The deflection profile at x=0 (Fig. 22) is in excellent agreement with the FEM solution obtained by ABAQUS (with 256 quadratic quadrilateral plate elements and 833 nodes). The relative error of maximum deflection is 0.29\mathrm{\%}.

The numerical calculation was also made by using n=81 and m=1. The relative error obtained is {\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}= −0.52%, reflecting that the number of terms n can be effectively reduced by increasing interior points m.

From the above analyses, it can be seen that the present approach is efficient in solving complex problems of plate bending, and good accuracy can be achieved with less calculation. Since this approach uses the whole-domain defined basis functions, which satisfy the governing differential equation of plate bending exactly, the deflection function can be approximated as the linear combination of basis functions, simply and directly. The unknown coefficients in the approximated deflection can only be determined with boundary conditions, and due to this, no domain integration is required. Moreover, the final deflection and internal forces are given in analytical form, allowing easy calculation.

In comparison with the traditional meshless methods, the main difference lies in the choice of basis function. The traditional meshless methods use radial basis functions [12,13], which are locally defined and do not satisfy the governing differential equation of plate bending. Because of this feature, a technique, e.g., moving least square method, is needed to construct the deflection function. The construction process is usually complex and time-consuming. When the principle of virtual work is applied to consider equilibrium conditions, domain integration becomes necessary, leading to a high cost of calculation.

6 Conclusions

A new meshless approach is proposed for the bending analysis of thin plates. Through definition of a set of points in plate domain, this work constructs a series of basis functions, each of which satisfies the governing differential equation of plate bending exactly. Plate deflection is then approximated as the linear combination of those functions, while the unknown coefficients in the approximation are determined through boundary conditions by using a collocation point method. With this approach, the solutions of both deflection and internal forces are given in analytical form.

This approach is further applied to analyze plates with arbitrary shapes and boundary conditions subjected to various loads. Numerical results show that the approach is simple and efficient in solving complex problems of plate bending, and good accuracy can be achieved with less calculation.

Theoretically, the approach could be extended to solve other problems, for example, plane elastic problems. This work is still in progress.

References

[1] Ugural A C. Plates and Shells: Theory and Analysis. 4th ed. Boca Raton: CRC Press, 2018

[2] Oñate E. Structural Analysis with the Finite Element Method. Linear Statics: vol. 2: Beams, Plates and Shells. Barcelona: International Center for Numerical Methods in Engineering (CIMNE), 2013

[3] Katsikadelis J T. The Boundary Element Method for Plate Analysis. London: Elsevier, 2014

[4] Karttunen A T, von Hertzen R, Reddy J N, Romanoff J. Exact elasticity-based finite element for circular plates. Computers & Structures, 2017, 182 : 219– 226

[5] Nguyen-Xuan H. A polygonal finite element method for plate analysis. Computers & Structures, 2017, 188 : 45– 62

[6] Karttunen A T, von Hertzen R, Reddy J N, Romanoff J. Shear deformable plate elements based on exact elasticity solution. Computers & Structures, 2018, 200 : 21– 31

[7] Katili I, Batoz J L, Maknun I J, Lardeur P. A comparative formulation of DKMQ, DSQ and MITC4 quadrilateral plate elements with new numerical results based on s-norm tests. Computers & Structures, 2018, 204 : 48– 64

[8] Videla J, Natarajan S, Bordas S P A. A new locking-free polygonal plate element for thin and thick plates based on Reissner-Mindlin plate theory and assumed shear strain fields. Computers & Structures, 2019, 220 : 32– 42

[9] Mishra B P, Barik M. NURBS-augmented finite element method for static analysis of arbitrary plates. Computers & Structures, 2020, 232 : 105869–

[10] Nhan N M, Nha T V, Thang B X, Trung N T. Static analysis of corrugated panels using homogenization models and a cell-based smoothed mindlin plate element (CS-MIN3). Frontiers of Structural and Civil Engineering, 2019, 13( 2): 251– 272

[11] Liu G R. Meshfree Methods: Moving Beyond the Finite Element Method. 2nd ed. Boca Raton: CRC Press, 2010

[12] Leitão V M A. A meshless method for Kirchhoff plate bending problems. International Journal for Numerical Methods in Engineering, 2001, 52( 10): 1107– 1130

[13] Liu Y, Hon Y C, Liew L M. A meshfree Hermite-type radial point interpolation method for Kirchhoff plate problems. International Journal for Numerical Methods in Engineering, 2006, 66( 7): 1153– 1178

[14] Chen J S, Hillman M, Chi S W. Meshfree methods: Progress made after 20 years. Journal of Engineering Mechanics, 2017, 143( 4): 04017001–

[15] Zhang H J, Wu J Z, Wang D D. Free vibration analysis of cracked thin plates by quasi-convex coupled isogeometric-meshfree method. Frontiers of Structural and Civil Engineering, 2015, 9( 4): 405– 419

[16] Pirrotta A, Bucher C. Innovative straight formulation for plate in bending. Computers & Structures, 2017, 180 : 117– 124

[17] Battaglia G, Di Matteo A, Micale G, Pirrotta A. Arbitrarily shaped plates analysis via Line Element-Less Method (LEM). Thin-walled Structures, 2018, 133 : 235– 248

[18] Nguyen V P, Anitescu C, Bordas S P A, Rabczuk T. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117 : 89– 116

[19] Zheng H, Liu Z J, Ge X R. Numerical manifold space of hermitian form and application to Kirchhoff’s thin plate problems. International Journal for Numerical Methods in Engineering, 2013, 95( 9): 721– 739

[20] Guo H W, Zheng H, Zhuang X Y. Numerical manifold method for vibration analysis of Kirchhoff’s plates of arbitrary geometry. Applied Mathematical Modelling, 2019, 66 : 695– 727

[21] Guo H W, Zhuang X Y, Rabczuk T. A deep collocation method for the bending analysis of Kirchhoff plate. CMC-Computers Materials & Continua, 2019, 59( 2): 433– 456

[22] Zhuang X Y, Guo H W, Alajlan N, Zhu H H, Rabczuk T. Deep autoencoder based energy method for the bending, vibration, and buckling analysis of Kirchhoff plates with transfer learning. European Journal of Mechanics-A/Solids, 2021, 87 : 104225–

[23] Guo H W, Zhuang X Y. The application of deep collocation method and deep energy method with a two-step optimizer in the bending analysis of Kirchhoff thin plate. Chinese Journal of Solid Mechanics, 2021, 42(3): 249– 266 (in Chinese)

[24] Li S C, Dong Z Z, Zhao H M. Natural Boundary Element Method for Elastic Thin Plates in Bending and Plane Problems. Beijing: Science Press, 2011 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press 2022.