Reduced-order strategy for meshless solution of plate bending problems with the generalized finite difference method

Abstract

This paper presents some recent advances on the numerical solution of the classical Germain-Lagrange equation for plate bending of thin elastic plates. A meshless strategy using the Generalized Finite Difference Method (GFDM) is proposed upon substitution of the original fourth-order differential equation by a system composed of two second-order partial differential equations. Mixed boundary conditions, variable nodal density and curved contours are some of the explored aspects. Simulations using very dense clouds and parallel processing scheme for efficient neighbor selection are also presented. Numerical experiments are performed for arbitrary plates and compared with analytical and Finite Element Method solutions. Finally, an overview of the procedure is presented, including a discussion of some future development.

Keywords
Plate bending; thin plates; generalized finite difference method; reduced-order; meshless; static

1 INTRODUCTION

The solution of differential equations by the Finite Difference Method (FDM) is a classic approach in regularly distributed nodes. It is a simple, intuitive and universal numerical method that deals directly with the differential form of the problem. However, disadvantages arise about the application of boundary conditions and treatment of arbitrary grids, required for the solution of irregular geometries. With the development of the Finite Element Method (FEM), versatile and able to easily deal with the problems mentioned above, scientific interest of researchers in new procedures associated to the FDM gradually diminished. Despite that, alternative versions of the classical form of this method were successfully developed in the 1960s to the 1990s. Advances arose in relation to the solution of arbitrary grids with the generalization of the method (Jensen, 1972 and Perrone and Kao, 1975), formulation of macroelements (Ghali and Bathe, 1970 and Bhattacharya, 1986), applications to nonlinear problems (Kao and Perrone, 1972), formulations based on energy principles (Pavlin and Perrone, 1979 and Fielding et al., 1997), among others .

In the last two decades, the possibility of unstructured point clouds with the generalized version of the finite difference method (GFDM) achieved an important status, with possible applications to adaptive refinements in h and p versions, with some examples obtained by Benito et al. (2002) and Gavete et al. (2018), who explored the reduction of local errors by adding new degrees of freedoms to the existing cloud. An h-adaptive method using a frame decomposition approach with conventional finite differences was also proposed by Srinivasa (2006).

Unlike the FEM, GFDM has no problems with hanging nodes and allows for cloud density transitions and the inclusion of new points by suitable operators. For this reason, it is also defined as the Meshless Finite Difference Method (MFDM). Thus, it is also natural the presence of recent works that seek the coupling between FEM and GFDM, combining the individual advantages of each method (Jaśkowiec and Milewski, 2016).

GFDM offers great flexibility in selection of a stencil for a given derivative. The final form of this selection is usually described as a star. On the other hand, random criteria may lead to singularity or bad conditioning of the coefficient matrix, limiting the quality of the derivatives. As described by Liu (2010), automation of nodal selection and improvement of the stability of the solution are still some challenges in this method, which does not require any kind of background mesh.

The works of Forsythe and Wasow (1960), Collatz (1966) and Jensen (1972) are considered pioneers in the solution of differential equations by finite differences in irregular grids. In all cases, five-point neighbors were used to approximate derivatives up to the second order. The GFDM method was successfully applied to large deflection of membranes (Kao and Perrone, 1972, Perrone and Kao, 1975), heat transfer problems (Liszka and Orkisz, 1980 and Chung, 1981), torsion problems and free vibration of membranes (Pulino Filho, 1989). For higher order derivative approximation, a larger number of neighbors is required, increasing the chances of a singular coefficient matrix. These difficulties were reported by Tseng and Gu (1989), who proposed a strategy of successive approximation of fourth order derivatives using order two operators, bypassing the problem of 14 neighbors in the Germain-Lagrange plate bending equation. Virtual points were produced to impose boundary conditions on square and elliptical plates.

Scientific papers describing GFDM's efforts in random point distributions are mostly focused on second-order partial differential equations or in a system composed of such equations. Node selection remains an active research field with the work of Benito et al. (2001) describing a study of several factors in the GFDM’s stars.

The singularity of the coefficient matrix has been bypassed by Prieto et al. (2011), who solved differential equations with 26-node stars and this development was also followed by Ureña et al. (2011), who solved bending of thin and thick plates, with the application of stars of up to 24 nodes.

Nowadays, some of the most recent advances include: computational acoustics (Wei et al., 2015), the solution of nonlinear water waves using potential flow theory (Zhang et al., 2016 and Fan et al., 2018), application to bidimensional shallow water equations (Li and Fan, 2017) and tridimensional adaptive cloud refinement (Gavete et al., 2018). As described by Zhang et al. (2016), the GFDM is versatile enough for many engineering applications. These recent contributions also demonstrate that the method is being rediscovered by researchers worldwide.

In an overview, there are very few papers dealing with plate bending combined with the GFDM, especially when describing the caveats and procedures for accurate results on plates with arbitrary geometries, including the computation of internal forces, which requires detailed investigation for irregular distributed nodes. As described previously, the stability of fourth-order operator is still a challenge in this method, and the vast majority of works dealing with this problem are limited to displacement evaluation, without further considerations on the precision of the field derivatives.

In a view of some of the above described problems the scope of the present work is given by five major contributions: (i) solution of the fourth order Germain-Lagrange partial differential equation, which is conveniently replaced by the simultaneous solution of two second-order equations (Marcus, 1932) avoiding higher order operators and singularities in the coefficient matrix; (ii) a technique for applying directional derivatives in curved contours; (iii) a parallel processing scheme for efficient neighbor selection in very dense clouds; (iv) a study of plate internal forces using uniform and varying density cloud nodes; (v) some recent results with automatic cloud generation including preprocessing of virtual nodes.

To the authors' knowledge, these are original and unexplored features in the GFDM applied to plate bending problems, which allows the solution of mixed boundaries in curved contours, in a fourth order differential equation, with the application of traditional preprocessors and the same simple background of the FDM, using a system of second-order partial differential equations. Details of the formulations and application examples are presented, followed by validations with theoretical and numerical values, demonstrating the highlights, advantages and disadvantages of the proposed procedures.

2 ORDER REDUCTION IN THE GERMAIN-LAGRANGE PLATE BENDING PROBLEM

The Germain-Lagrange equation in its classical form applied to plates of uniform stiffness and small thickness is given by:

$\frac{{\partial }^{4}w}{\partial x^4}+2\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+\frac{{\partial }^{4}w}{\partial y^4}=\frac{q(x,y)}{D}$ (1)

which may have its order reduced by a system of second-order partial differential equations (Timoshenko and Woinowsky-Krieger, 1959):

$\{\begin{array}{c}\frac{{\partial }^{2}M}{\partial x^2}+\frac{{\partial }^{2}M}{\partial y^2}=-q\left(x,y\right)\\ \frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}=-\frac{M}{D}\end{array}$ (2a,b)

where $w$,$q$ and $D$ denote the transverse displacement, the applied load and the flexural rigidity of the plate, respectively. Finally, $M$corresponds to an equivalent bending moment given by:

$M=\frac{M_x+M_y}{1+v}=-D\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right)$ (3)

with $M_x$ and $M_y$defining the bending moments in directions $x$ and $y$, respectively. The term $v$denotes the Poisson´s ratio. At the boundary $\Gamma$ it follows:

$M_{\Gamma }=-D\left(\frac{{\partial }^{2}w}{\partial n^2}+\frac{{\partial }^{2}w}{\partial t^2}\right)$ (4)

with $\overrightarrow{t}$and $\overrightarrow{n}$ defining the tangent and normal directions to the contour, as shown in Fig. 1, and α is the angle between the horizontal direction and $\overrightarrow{n}$.

Figure 1
Normal and tangential directions at the contour.

The advantage of the order reduction approach is evident since an originally fourth-order problem is reduced to the simultaneous solution of two Poisson equations (which is a widely studied equation in GDFM related works). In the scenario of a simply supported polygonal plate, an additional advantage arises with the successive solution of this system. Therefore, one can solve (2a) for internal values of $M$since at the boundary:

$M_{\Gamma }=-D(\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^{2}w}{\partial n^2})=0$ (5)

$w_{\Gamma }=0$ (6)

and latter substitute these values in (2b) for the displacement solution. This probably provides one of the simplest and most elegant solution for the Germain-Lagrange problem and has been extensively described and applied by many authors, including Bhattacharya (1986) and Ugural (1981). It is observed, however, that condition (5) is not satisfied for curved contours (Timoshenko and Woinowsky-Krieger, 1959) and boundary moments are unknown at these scenarios. Therefore, the successive strategy is restricted to simply supported polygonal plates. In a general scenario of plates with curved contours, the system must be solved simultaneously.

For most general cases, contour conditions for clamped (Guang-yao and Han-bin, 1993) and simply supported boundaries (Tseng and Gu, 1987) are given respectively by:

${\left[\frac{\partial w}{\partial n}\right]}_{\Gamma }=\frac{\partial w}{\partial x}\cos \left(\alpha \right)+\frac{\partial w}{\partial y}sin\left(\alpha \right)=0$ (7)

${\left[\frac{{\partial }^{2}w}{\partial n^2}\right]}_{\Gamma }=\left(\frac{{\partial }^{2}w}{\partial x^2}+v\frac{{\partial }^{2}w}{\partial y^2}\right){\cos }^2\alpha +\left(\frac{{\partial }^{2}w}{\partial y^2}+v\frac{{\partial }^{2}w}{\partial x^2}\right){\sin }^2\alpha +2(1-v)\frac{{\partial }^{2}w}{\partial x\partial y}\sin \alpha \cos \alpha =0$ (8)

The simultaneous solution requires a larger number of unknowns since displacements $w$and bending moments$M$are solved simultaneously. Still, the procedure is straightforward with finite difference operators applied to (2a) and (2b). Once solution for displacements is achieved, bending moments are evaluated by displacement field derivatives:

$M_x=-D\left(\frac{{\partial }^{2}w}{\partial x^2}+v\frac{{\partial }^{2}w}{\partial y^2}\right)$ (9)

$M_y=-D\left(\frac{{\partial }^{2}w}{\partial y^2}+v\frac{{\partial }^{2}w}{\partial x^2}\right)$ (10)

$M_{xy}=D(1-v)\left(\frac{{\partial }^{2}w}{\partial xy}\right)$ (11)

3. THE BASIC IDEA OF GFDM AND SOME PROPOSED PROCEDURES

GFDM maintains the simplicity of partial derivative approximations as in FDM and at the same time uses points in the domain without prior relationship between them. This aspect allows great freedom in the choice of differential operators and treatment of arbitrary geometries. According to Liu (2010), two significant advantages are inherent to GFDM over other meshless category methods: ease of implementation and absence of integral formulations. The basic idea of this method arises with the fact that any differentiable function can be expanded in a bidimensional Taylor series. Consider the two-dimensional expansion of $f(x,y)$, centered in $P(x_0,y_0)$ at the domain $\Omega$, with $\Omega \subset {ℝ}^2$and boundary $\Gamma$, as depicted in Fig. 2. A general expansion is given by:

$f\left(x,y\right)=f_0+{\displaystyle \sum _{z=0}^{\infty }\left\{\frac{1}{z!}{\displaystyle \sum _{s=0}^z\left(\begin{array}{c}z\\ s\end{array}\right)}\frac{{\partial }^z{f}_o}{\partial x^{z-s}\partial y^{\text{}s}}{\left(x-x_o\right)}^{z-s}{\left(y-y_o\right)}^s\right\}}$ (12)

with

$\left(\begin{array}{c}z\\ s\end{array}\right)=\frac{z!}{(z-s)!s!};f_0=f\left(x_0,y_0\right)$ (13a,b)

and boundary conditions:

$m(f)+r=0\begin{array}{cc}& in\begin{array}{cc}& \Gamma \end{array}\end{array}$ (14)

with $m$ defining a linear differential operator, and$r$describing an arbitrary function.

Figure 2
(a) General scheme with boundary, contour and neighbors; (b)corresponding star

As presented by Perrone and Kao (1975) and many others, the evaluation of a given derivative requires a number of neighboring nodes equal in number to the total terms developed up to the desired order. Second-order approaches require at least five points, whereas fourth-order approximations use 14 neighbors. In general, a system of equations using (12) can be converted to matrix form with Eq. (15).

$\left[\begin{array}{cccccccc}{h}_1& k{}_1& \frac{h_1{}^2}{2}& h_1{k}_1& \frac{k_1{}^2}{2}& \frac{h_1{}^3}{6}& \frac{h_1{}^2{k}_1}{2}& \cdots \\ h_2& k{}_2& \frac{h_2{}^2}{2}& h_2{k}_2& \frac{k_2{}^2}{2}& \frac{h_2{}^3}{6}& \frac{h_2{}^2{k}_2}{2}& \cdots \\ h_3& k{}_3& \frac{h_3{}^2}{2}& h_3{k}_3& \frac{k_3{}^2}{2}& \frac{h_3{}^3}{6}& \frac{h_3{}^2{k}_3}{2}& \cdots \\ h_4& k{}_4& \frac{h_4{}^2}{2}& h_4{k}_4& \frac{k_4{}^2}{2}& \frac{h_4{}^3}{6}& \frac{h_4{}^2{k}_4}{2}& \cdots \\ h_5& k{}_5& \frac{h_5{}^2}{2}& h_5{k}_5& \frac{k_5{}^2}{2}& \frac{h_5{}^3}{6}& \frac{h_5{}^2{k}_5}{2}& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\left\{\begin{array}{c}\frac{\partial f_0}{\partial x}\\ \frac{\partial f_0}{\partial y}\\ \frac{{\partial }^2{f}_0}{\partial x^2}\\ \frac{{\partial }^2{f}_0}{\partial x\partial y}\\ \frac{{\partial }^2{f}_0}{\partial y^2}\\ ⋮\end{array}\right\}=\left\{\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\\ f_5\\ ⋮\end{array}\right\}-\left\{\begin{array}{c}{f}_0\\ f_0\\ f_0\\ f_0\\ f_0\\ ⋮\end{array}\right\}$ (15)

With $i$ denoting a local neighbor, $h_i=x_i-x_0$ and$k_i=y_i-y_0$. Eq. (15) can be further represented by a compact matrix notation using an approximation with $N$ neighbors:

$A_{{}_{N\times N}}{u}_{{}_{N\times 1}}=f_{{}_{N\times 1}}-f_0{}_{{}_{N\times 1}}$ (16)

The partial derivatives in $u$ are obtained with:

$u=Bf-B{f}_0;B=A^{-1}$ (17a,b)

The generalized form of any derivative is given by:

$u_i={\displaystyle \sum _{j=1}^N{b}_{ij}\left(f_j-f_0\right)}$ (18)

where $u_i$and $b_{ij}$are indexes of the following matrices:

$u={\left[\begin{array}{ccccc}\frac{\partial f_0}{\partial x}& \frac{\partial f_0}{\partial y}& \frac{{\partial }^2{f}_0}{\partial x^2}& \frac{{\partial }^2{f}_0}{\partial x\partial y}& \cdots \end{array}\right]}^T={\left[\begin{array}{ccccc}{u}_1& u_2& u_3& u_4& \cdots \end{array}\right]}^T;B=\left(b_{ij}\right)\in {ℝ}^{N\times N}$ (19a,b)

3.1. Selection of neighboring nodes

The choice of neighboring nodes in the GFDM can be based on several criteria. The simplest option (and that should be discouraged) is the selection of the nodes closest to P (distance criterion). As shown by Perrone and Kao (1975), problems can arise due to point alignment, producing singularity in A. In view of the previous problem, some rules were proposed for adequate selection of the five neighbors. One of them consists on a quadrant division with origin at P (Liszka and Orkisz, 1980 and Benito et al., 2001), where each subregion must be occupied by a node (Fig 3a). Another, proposed by Perrone and Kao (1975), is given by the arrangement in circular pie-type sectors (Fig 3b). In both cases, the goal is the equidistribution of points in a star-like configuration, which arises because of the geometric appearance formed around a central node (as shown in Fig. 2b). In a way, the problem consists on the equilibrium of two important variables: good distribution around the central node P and smaller possible distances. It is also important to emphasize that very strict criteria should be avoided, as they could: (i) fail on capture neighboring points; or (ii) become excessive time demanding on dense clouds. The algorithms developed in the present work use a process of three simultaneous selection criteria by means of (Fig. 3c): Quadrants (Q), Distances (d); Minimum angles ($\theta$).

Figure 3
(a) four quadrant criterion, (b) pie sector criterion, (c) actual proposal;

The last criterion is a strategy where the angles between the lines connecting the central node to the neighboring points are computed. An angular tolerance is established and avoids points aligned or very close in the same quadrant or at the border of two adjacent quadrants. These steps are presented in Table 1. The authors further indicate that expanded versions of the procedure can be obtained for a greater number of points. For stars with nine neighbors, for example, one could just duplicate the total of quadrants available from Step 1.

Table 1
Proposed strategy for selection of 05 internal neighbors.

Step	Procedure
1	Four quadrants are organized as Q1-2-3-4 and an angular tolerance is prescribed.
2	Points are sorted in ascending order with respect to distance to P. A flag with the corresponding quadrant is also prescribed.
3	Selection of the closest point to P. The corresponding quadrant is removed from Step 1
4	The second point is the next based on the available quadrants and under the prescribed angular tolerance. The corresponding quadrant is removed from Step 1.
5	The procedure on Step 4 is repeated until the selection of the fourth point. The vector containing the list of available quadrants is now empty.
6	Step 2 is restarted without the four selected nodes. The fifth and last point is the closest to P that guarantees the angular tolerance among the other selected neighbors. No restriction is applied for the corresponding quadrant.

4. COMPUTATIONAL ASPECTS

The routines of the present work were developed with MATLAB 2017a software. The code is structured for the solution of arbitrary plates with curved contours and mixed boundary conditions (simply supported and clamped). The main aspects are presented below.

4.1. Cloud generation

Point clouds are constructed from two distinct paths:

1. Proprietary generator, which bombards the domain interior with random points and where a tolerance (in the distance) is established between the next point to be inserted and those already existing (Fig. 4a). The procedure is repeated according to user-defined iterations. Regions with different nodal densities are also possible;

2. GiD v.11.0.7. using Delaunay triangulation. In both cases, only the nodal coordinates are recovered. Strategy 2 is extremely useful for complex computational domains (containing holes, irregular) and with diverse cloud density transitions since generator described in 1 is still in the early stages of development. An additional routine is applied to construct virtual points in both strategies. The procedure is relatively simple and associates a virtual point with each point in the boundary. Initially, the mean normal direction of the vertices adjacent to a point in the contour is computed, and a virtual point is projected with arbitrary distance δ in this direction (usually preserving the same order of magnitude of the contour discretization to its nearest neighbors within the domain). Fig. 4b illustrates this procedure. The algorithm developed by Kroon (2011) was applied at this preprocessing stage, which calculates curvatures and normal directions of lines and contours. In this case, quadratic functions are employed in a polynomial fit of the contour and finite differences are applied for calculation of the normal direction.

Figure 4
(a) Automatic cloud generation and (b) virtual nodes with distance δ.

4.2. Neighbor selection with parallel processing

An interesting aspect of the current proposal is the use of parallel processing for neighbor selection. For very dense clouds the time spent in building a matrix with the best neighbors at each point of the domain becomes a tedious operation. The classification of distances about the point P is made with MATLAB’s QuickSort algorithm, which presents O (n log n) complexity. For domains with a large number of 'n' points, this step is a great disadvantage of GFDM since no triangularization or prior ordering between the points is defined. On the other hand, the use of parallel processing appears as an interesting and accessible workaround on multi-core computers. In MATLAB this procedure only requires the use of a simple PARFOR command, which performs simultaneous operations on the number of processors allocated to a task. In the present work desktops and laptops with four physical cores were used. In this case, the domain is partitioned into four subregions (no significant advantage is observed with the use of hyperthreading).

In quantitative terms, Fig. 5 presents the processing time and the performance gain obtained by a parallel processing option. The term speedup refers to the ratio between the processing time of the Single (1 core) and Quad (4 cores) options. The simulations of this example were performed in a square domain with an Intel i7 7700HQ processor (four cores, eight threads, maximum frequency 3.8GHz). Expressive results are observed and can be improved with the use of more robust computers or distributed processing.

Figure 5
Performance output: (a)clouds, options and processing time; (b) speedup using parallel processing.

An alternative to a parallel processing option arises through the p-adaptive refinement of the star. However, this requires even more efficient neighbor selection routines. For the moment the authors discourage this type of alternative due to the possibility of sigularity in matrix A. For convenience the neighbors are stored in indexes of the following matrix:

$V={\left[V_1\mathrm{,...,}{V}_l\mathrm{,..,}{V}_S\right]}^T\in {ℝ}^{S\times 5};V_l=\left[C_1\mathrm{,...,}{C}_5\right]$ (20a,b)

where $V_l$denotes the corresponding line $l$for the five neighbors C´s of a point P.

4.3. Generalized finite difference operators

The used operators arise from Eq. (18). In this work, only five neighbors are used for each point P, therefore $j=\mathrm{1..5}$. For a second-order derivative, e.g. $i=3$:

$u_3\approx {\displaystyle \sum _{j=1}^5{b}_{3j}}{f}_j-{\displaystyle \sum _{j=1}^5{b}_{3j}}{f}_o$ (21)

The system of linear equations is assembled after storing the neighbor’s information. Here the routine is divided in two directions: (i) successive or (ii) simultaneous solution. Procedures are presented in the next section.

4.4. Matrix assembly procedures - successive solution

This strategy stands out for its simpler implementation and lower computational cost when compared to the simultaneous solution strategy. Using the second order operators in Eq. (2a), and adopting $f\to M$:

$\left(b_{3j}+b_{5j}\right)M_I-M_l{\displaystyle \sum _{j=1}^5\left(b_{3j}+b_{5j}\right)}=-q_l,\begin{array}{cc}& l=\mathrm{1..}ni\end{array};$ (22)

$K{1}_{ni\times ni}{M}_{ni\times 1}=-q_{ni\times 1}$ (23)

Where $l$ denotes the corresponding point P and $ni$ the number of internal points of the domain. It should be noted that $M_I$ is the vector with the equivalent moments in the global numbering of the neighboring points, being:

$M_I{}^T={\left\{\begin{array}{ccc}{M}_{C_1},& \mathrm{...}& ,M_{C_5}\end{array}\right\}}^T$ (24)

Once the solution for $M$is achieved and following the same previous procedure, now using Eq. (2b) and with $f\to w$:

$\left(b_{3j}+b_{5j}\right)w_I{}_{{}_{5\times 1}}-w_l{\displaystyle \sum _{j=1}^5\left(b_{3j}+b_{5j}\right)}=-\frac{1}{D}{M}_l,\begin{array}{cc}& l=\mathrm{1..}ni\end{array};$ (25)

$K{2}_{ni\times ni}{w}_{ni\times 1}=-\frac{1}{D}{M}_{ni\times 1}$ (26)

Where $w_I$ is the vector with the displacements in the global numbering of the neighboring points, being:

$w_I{}^T={\left\{\begin{array}{ccc}{w}_{C_1},& \mathrm{...}& ,w_{C_5}\end{array}\right\}}^T$ (27)

Matrices $K1$ and $K2$ are nonsymmetric and stored using sparse capabilities. Vectors $M$and $w$include all equivalent moments and displacements.

4.5. Matrix assembly procedures - simultaneous solution

The values of $M$are not initially known in the contour ($M_{\Gamma }\ne 0$in contrast to the previous case). Thus, an option is made for a simultaneous solution of the system of equations. The first step consists on application of (22). With this$ni$equations are generated, where$ni$corresponds to the total of internal points.

At the end of the first step, there are unknowns in number $npt=ni+nc$ (total number of points, internal and contour), since the values of $M$at the contour are additional unknowns. The next step is given by application of:

$\left(b_{3j}+b_{5j}\right)w_I{}_{{}_{5\times 1}}-w_l{\displaystyle \sum _{j=1}^5\left(b_{3j}+b_{5j}\right)}+\frac{M_l}{D}=\mathrm{0,}l=\mathrm{1..}npt;$ (28)

or in a compact notation:

$K{2}_{npt\times (ni+nv)}{w}_{(ni+nv)\times 1}+\frac{1}{D}{M}_{npt\times 1}=0$ (29)

Equation (29) is valid for real points on the domain, which implies in unknowns of displacements of the virtual points. For the bending moments, the same does not occur. Table 2 presents a summary of the involved parameters.

Table 2
Summary of different strategies on the simultaneous solution.

Scenario	P range for $l$	Unknowns in $M$	Unknowns in $w$
Simultaneous solution	1..ni for $M$ 1..npt for $w$	Real nodes	Real and virtual nodes

Thus $nc$new equations are obtained. Here the importance of the preprocessing step with $nc=nv$ is verified (Fig. 4b), where $nv$denotes the number of virtual nodes. Otherwise, the system will be undetermined. The next section describes the required procedures.

4.5.1 Boundary conditions

For convenience, only cases with simply supported boundary, clamped or a combination of both (mixed) are addressed. Thus, the condition in (6) is always satisfied.

Application of differential operators in (7) for the clamped contour results in:

$\begin{array}{l}{b}_{1j}{w}_I\cos \alpha +b_{2j}{w}_I\sin \alpha =\mathrm{0,}\begin{array}{cc}& l=ni+\mathrm{1..}npt\end{array};\\ \end{array}$ (30)

and following from (8) in the case of a simply supported boundary:

$\begin{array}{l}\left(b_{3j}{w}_I+v{b}_{5j}{w}_I\right){\cos }^2\alpha +\left(b_{5j}{w}_I+v{b}_{3j}{w}_I\right){\sin }^2\alpha +\\ 2(1-v)b_{4j}{w}_I\sin \alpha \cos \alpha =0\begin{array}{cc}& l=ni+\mathrm{1..}npt;\end{array}\end{array}$ (31)

For a combination of boundary conditions (simply supported and clamped), the range $l$ should be modified according to the problem. In a compact notation:

$K{3}_{nc\times (ni+nv)}{w}_{(ni+nv)\times 1}=0$ (32)

4.6. Resulting linear system

In the case of a successive solution, the procedure is immediate and consists of the solution from (23) to $M$and after (26) for $w$. For the simultaneous solution the following linear system is defined:

$\left[\begin{array}{cc}K{1}_{ni\times npt}& \\ {\overline{M}}_{ni\times ni}& K{2}_{npt\times (ni+nv)}\\ & K{3}_{nc\times (ni+nv)}\end{array}\right]X=F;\overline{M}=\frac{1}{D}\left[\begin{array}{cccc}{M}_1& & & 0\\ & M_2& & \\ & & \ddots & \\ 0& & & M_{ni}\end{array}\right]$ (33a,b)

$X={\left(\begin{array}{cccccccccccccccccc}{M}_1& M_2& \dots & M_l& \dots & M_{npt}& w_1& w_2& \dots & w_l& \dots & w_{ni}& {\widehat{w}}_1& {\widehat{w}}_2& \dots & {\widehat{w}}_l& \dots & {\widehat{w}}_{nc}\end{array}\right)}^T$ (34)

$F=-{\left(\begin{array}{ccccccccc}{q}_1& q_2& \dots & q_l& \dots & q_{npt}& 0& 0& \dots \end{array}\right)}^T$ (35)

with $\widehat{w}$ denoting the displacements of the virtual nodes.

5.0 Numerical results

5.1. Case 01 - Simply supported square plate with downward point load

Since this is a simply supported plate, the successive solution of two differential equations presents an elegant strategy and reduces the problem to second-order equations with first-type boundary conditions, as described on the computational procedures. The following example (Fig. 6a) is given by a simply supported square plate with an elastic modulus $E=20GPa$, Poisson coefficient $v=0.2$ and thickness $h=0.05m$. A downward point load $P={10}^{6}N$ is applied at the plate center. Since the GFDM requires distributed loading, the point load is divided by ​​the influence region of the central point. This influence domain is evaluated using Delaunay triangulation of the points near this region. The central node receives a third of the areas of all the triangles to which it is connected, and the equivalent distributed load is computed with the division of $P$by this area.

Figure 6
Plate geometry and variable density point clouds for Case 01 (RV-type).

As described by Szilard (2004) the following closed-form solution is available:

$w(x,y)=\frac{4P}{{\pi }^4{L}^{2}D}{\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }\frac{\sin (m\pi /2)\sin (n\pi /2)}{{\left[\left(\frac{m^2}{L^2}\right)+\left(\frac{n^2}{L^2}\right)\right]}^2}}}\sin \frac{m\pi x}{L}\sin \frac{n\pi y}{L},\begin{array}{cc}& m,n=\mathrm{1,2,3}\dots \end{array}$ (36)

where $L$ denotes the plate dimension. Internal forces are evaluated with application of (9)-(11). Table 3 presents a summary of the analyzed clouds. The uniform clouds form an orthogonal grid, whereas clouds with varying densities are shown in Figs. 6b-d.

Table 3
Summary of analyzed clouds and properties (Case 01).

Cloud	Cloud node density type	Number of nodes	Boundary condition
RU1	Uniform (RU)	441	Simply supported
RU2		1681
RU3		2601
RV1	Variable (RV)	585
RV2		1887
RV3		2871

Results of displacements in section cut R are shown in Fig. 7. Convergence of this parameter is possible at the center of the plate (albeit slowly) and occurs more quickly for the RU type clouds. On the other hand, the convergence of internal forces is impossible even for the closed-form solution at the center of the plate (Szilard, 2004). This will be highlighted in the next comparisons.

Figures 8a-c show the behavior of the uniform cloud set RU with respect to the bending moment Mx at section cut R. In this case the exact solution was computed with 10⁴, 10⁶ and 10⁸ terms in the double series. It is worth noticing the rapid convergence of the results at points that are in about 90% of the domain. In spite of that, convergence is impossible at the plate center. The exact solution displays a singularity at this point.

In Figs. 8d-f clouds with variable node densities RV are analyzed. In a similar fashion to RU clouds, results are also in excellent agreement with the closed-form solution. The localized refinement at the plate center improves results at the center of the domain, and this becomes clear when Figs. 8b and 8e are compared. This same behavior is also noticeable for Figs. 8c and 8f. RV clouds showcase the great versatility of GFDM when dealing with singularities since no specific transition is required.

Figure 7
Displacement results for RU-type and RV-type clouds at section cut R.

Figure 8
Bending moment (Mx) results for RU-type and RV-type clouds at section cut R.

5.2. Case 02 - Simply supported octagonal plate with uniform load

Reference axis and analyzed sections are presented in Fig. 9a. Some examples of clouds with uniform point density are shown in Figs. 9b-d. Analysis parameters are given by an elastic modulus $E=20GPa$, Poisson coefficient $v=0.2$, thickness$h=0.25m$ and transverse load $q=1000Pa$.

Table 4 indicates all analyzed cases for this geometry. The PV-type cloud mimics an adaptive refinement (although no special criterion is employed) and aims to evaluate the performance of the GFDM for point density transitions. Fig. 10 presents these clouds.

Table 4
Summary of analyzed clouds and properties (Case 02).

Cloud	Cloud node density type	Number of nodes	Boundary condition
PU1	Uniform	234	Simply supported
PU2		422
PU3		1134
PU4		1586
PU5		2441
PU6		4406
PV1	Variable	1109
PV2		1501
PV3		2427

Figure 9
(a) Plate geometry and (b) uniform point clouds for Case 02 (PU-type).

Figure 10
Variable density point clouds for Case 02 (PV-type).

5.2.1 Results for PU-type clouds

The reference solution arises with a FEM model using quadrilateral plate bending elements with the same number of nodes as the PU6 cloud. Fig. 11 presents a series of vertical displacements along this plate for sections S1 and S4.

An additional investigation arises with the evaluation of Mx bending moments along the selected sections, as shown in Fig. 12. These results are excellent, showing discontinuities only at the polygon’s vertices (present at the extremities of sections S1 and S3). In these regions, the GFDM points a limitation in terms of the hypothesis of $M=0$, adopted in the sequential solution. In general, this value will be nonzero at the vertices and will be an additional unknown of the problem, even for simply supported plates. Nevertheless, the results are satisfactory for a polygon with 8 vertices.

Figure 11
Displacement results for PU-type clouds at S1 and S4 sections.

Figure 12
Bending moment results for PU-type clouds at S1 and S2 sections.

In global terms, the bending moments isovalues for the GFDM and FEM are presented in Fig. 13. Excellent agreement is observed.

Figure 13
Isovalues of bending moments Mx [Nm/m] for (a) GDFM and (b) FEM - Cloud PU6.

5.2.2. Results for PV-type clouds

Results for PV-type clouds are presented only for the bending moments isovalues, since an excellent agreement with FEM results is also observed for the vertical displacement field. Figure 14a illustrates these results. At first glance, these nonuniform clouds present a major drawback on the GFDM, since badly disposed stars (that are a good approximation for displacements) fail to capture the higher order derivatives involved on the bending moment computation, as shown in Fig. 14b (detail). However, this is majorly due to the bad arrangement provided by the PV1 cloud. As depicted in Fig. 15, some very nice results are achieved for clouds PV2 and PV3, which include both an aggressive and more conservative refinement towards the center of the polygon. In both scenarios, the boundary node density is the same. While PV1 refinement was achieved using the authors’ proposed node generator (which randomly bombards the domain, with no further restrictions except node distance) both PV2 and PV3 clouds where generated using GiD. Therefore, the problem is indeed involving and may require further restrictions, especially if accurate bending moments are required.

Figure 14
Isovalues of (a) displacements [m] and (b) bending moments Mx [Nm/m] - Cloud PV1.

Figure 15
Isovalues of bending moments Mx [Nm/m] - Clouds (a) PV2 and (b) PV3.

5.3. Case 03 - clamped elliptical plate with uniform load

In this example a clamped elliptical plate is analyzed, as depicted in Fig. 16a. Analysis parameters are given by an elastic modulus $E=20GPa$, Poisson coefficient $v=0.2$, thickness$h=0.05m$ and transverse load $q=1000Pa$. Figures 16b-d present some sample clouds and Table 5 summarizes the performed simulations.

Table 5
Summary of analyzed clouds and properties (Case 03).

Cloud	Cloud node density type	Number of nodes	Boundary condition
EU1	Uniform	432	Clamped
EU2		1005
EU3		3823
EU4		14901

Figure 16
(a) Plate geometry and (b,c,d) uniform point clouds for Case 03 (EU-type)

5.3.1 Results for EU-type clouds

In this case the simultaneous strategy described by Eq. (33a) must be applied. The reference solution is given by the resulting derivatives of the exact solution for displacements. It follows (Timoshenko and Woinowsky-Krieger, 1959):

$w_{exact}(x,y)=\frac{q}{D\left(\frac{24}{a^4}+\frac{24}{b^4}+\frac{16}{a^2{b}^2}\right)}{\left[1-{\left(\frac{x}{a}\right)}^2-{\left(\frac{y}{b}\right)}^2\right]}^2$ (37)

where $a$and $b$denote the ellipse’s outer and inner radius, respectively.

Figure 17a presents the bending moment along section cut E1, while the twisting moment is evaluated at section cut E2, as shown in Fig. 17c. A cumulative histogram for relative errors is provided for clouds EU1-EU3. In Fig. 17b the bending moment error is inferior to 10% for 56%, 82% and 92% of the total number of nodes of the analyzed clouds. As expected, higher node density provides better results. This same analysis is performed for the twisting moment (in Fig. 17d), resulting in values of 73%, 85% and 93% of the total nodes with errors inferior to 10%. This excellent agreement for clouds with less than 4,000 points indicates that the GFDM scales well on modern computers, with neighbor finding routines performing in less than ten seconds.

Figure 17
Section cuts and cumulative histograms for relative errors of $M_x$and$M_{xy}$ on the domain.

Figures 18a and 18b specify the percentage error, concerning the exact solution, for the GFDM and the FEM (with triangular plate bending elements) for the bending moments Mx and My, respectively. A slight advantage for the GFDM is observed for convergence of Mx at the center of the plate (Fig. 18a). On the other hand, for the bending moment My the FEM demonstrates a practical advantage (Fig. 18b). In both methods, the errors in bending moments are less than 0.16% for the EU3 cloud. For the twisting moment Mxy (Fig. 18c), the advantage of the FEM is remarkable, yet the errors in the GFDM are acceptable, reaching approximately 2% in the EU3 cloud.

Figure 18
Convergence study for bending moments at different point locations

5.4. Case 04 - irregular plate with mixed boundary (clamped and simply supported)

The following example is characterized by an irregular plate with mixed boundary conditions, as shown in Fig. 19a. An example of preprocessing for IU2 cloud with 913 points is also presented in Fig. 19b, including the required virtual nodes. Analysis parameters are given by an elastic modulus $E=21.3GPa$, Poisson coefficient $v=0.2$, thickness$h=0.80m$ and transverse load $q=100kPa$. Table 6 summarizes the performed simulations.

Table 6
Summary of analyzed clouds and properties (Case 04).

Cloud	Cloud node density type	Number of nodes	Boundary condition
IU1	Uniform	314	Mixed
IU2		913
IU3		1781
IU4		3336
IU5		9408
IU6		12750

Figure 19
(a) Irregular plate with mixed boundary conditions and (b) preprocessing example for IU2 cloud.

The reference solution is given by a FEM model using quadrilateral plate bending elements with the same number of nodes as the IU6 cloud. For the GFDM model a simultaneous solution scheme is applied since this example includes a clamped portion of the boundary. Here the author’s point an important consideration: in general, polygonal approximations will result in errors at the polygon’s vertices (as discussed on Case 1). Therefore, even if the current example was given by an entirely simply supported boundary, the successive strategy would still be insufficient. This irregular plate geometry is based on an ‘n’ sided polygon and the equivalent bending moments $M$ at these vertices are nonzero. Poor results are expected whenever ‘n’ increases, since these moments where assumed as zero. On the other hand, very nice quality approximations are expected whenever straight sides are present. The successive strategy overcomes these limitations.

5.4.1 Results for IU-type clouds

Figure 20a presents a convergence analysis of peak displacements using GFDM and the FEM reference results. The factor ‘rd’ denotes the relative difference between both methods. GFDM converges slowly when compared to FEM, but these differences are tolerable, being under 3% for IU3 cloud. It should be noticed that the simply supported boundary requires second-order and mixed derivatives, as shown in Eq. (8), and these will converge slowly when compared to problems with clamped boundaries.

Clouds IU1-IU5 are selected for displacement evaluation and IU2-IU6 are applied for computation of internal forces. The prescribed section cut is shown in Fig. 19a and these results are depicted in Figs. 20b-20d. Convergence at the simply supported boundary (pointed in Fig. 20c) is more sensible to point refinement and requires at least IU6 cloud for proper representation of the bending moments. In general, this is counteracted with dense clouds using the proposed parallel processing procedures, which greatly reduces the computational time.

Figure 20
(a) Convergence study for peak displacements, (b), (c) and (d) section cut results.

In global terms, the bending moments isovalues for the GFDM and FEM are presented in Figs. 21-22 for IU6 cloud. Excellent agreement is observed.

Figure 21
Isovalues of bending moments Mx [kNm/m] for (a) FEM and (b) GFDM - Cloud IU6.

Figure 22
Isovalues of bending moments My [kNm/m] for (a) FEM and (b) GFDM - Cloud IU6.

A detailed MATLAB code for these simulations is provided in the following GitHub link: https://github.com/GFDM-plate/scripts

6.0 Concluding remarks

This work presents a systematic solution for the bending of thin elastic plates based on the Generalized Finite Difference Method (GFDM). In this method, the mesh entity is non-existent (even the background mesh). Points can be moved, included, or removed with great ease in any region of the domain, without requiring the reconstruction of meshes or transitions between them. A system of two second-order equations replaces the original fourth-order partial differential equation. This allows a considerable advantage on simply supported plates with straight edges, which are solved in a successive strategy. It also enables a simple star configuration, with only five neighbors, reducing the chances of a singular coefficient matrix. A convenient procedure for neighbor selection using triple criteria: distance, angles, and quadrants is discussed and implemented using a parallel processing strategy with multi-core computers. This enables the solution of very dense point clouds, with up to 50,000 nodes. To the authors’ knowledge, this reduced-order strategy in simply supported plates with the GFDM is perhaps one of the most elegant techniques for polygonal domains, based on the same simple background of the conventional Finite Difference Method (FDM), extending the application range for arbitrary grids. However, there are some drawbacks. The successive strategy fails for multiple sided polygonal domains, since equivalent moments are nonzero at these vertices. On these scenarios, the successive strategy is not recommended. GFDM also relies on stability based on cloud arrangement and field derivatives are sensible to node density transitions. This is hardly a problem in displacements computations, but it is an issue whenever internal forces are required. Further studies are recommended in these cases.

Examples with curved boundaries and multiple boundary conditions were discussed. For these scenarios an additional strategy is proposed, where the displacements and equivalent moments are solved simultaneously. A technique for the inclusion of virtual nodes in conventional preprocessors is presented, based on the projection of boundary nodes at the normal contour direction. Directional derivatives were implemented for the application of boundary conditions on clamped and simply supported boundaries. Results for elliptical and irregular plates are in excellent agreement with reference results, and these studies were also extended for the analysis of internal forces. Finally, the analyzed examples show that the proposed algorithm maintains the simplicity inherent to the conventional FDM and enables the geometric generality for arbitrary plates, a characteristic found in the FEM. A good basis is set forth for the solution of other differential equations and additional plate bending problems with this method, which can be expanded to free vibration, nonlinear analysis, and many others.

In the scenario of layered plates with distinct materials (laminates) the procedure is relatively straightforward since an equivalent flexural stiffness (Szilard, 2004) is established as prescribed by Pister and Dong (1959). In this way, the order reduction procedure (proposed in the manuscript) remains valid. However, the strategy provided on this paper fails for orthotropic plates. As provided by Huber (1923) and Szilard (2004) a fourth order differential operator might be required, and this might lead to a singularity in the coefficient matrix. Therefore, for this specific case the authors recommend furthers studies in order to avoid uniqueness of matrix B.

In general, the FEM will surpass the GFDM in the most diverse applications. However, the extreme ease of formulation and the absence of any type of mesh make this numerical method extremely attractive for research beginners in meshless methods. In the authors' opinion, the combination of these two methods is promising and might result in a fascinating hybrid method in future publications (where loads and boundary conditions are applied to the FEM).

References

• Benito, J. J., Urena, F., Gavete, L., (2001). Influence of several factors in the generalized finite difference method. Applied Mathematical Modelling 25:1039-1053.
• Benito, J. J., Urena, F., Gavete, L., (2002). An h-adaptive method in the generalized finite differences. Computer Methods in Applied Mechanics and Engineering 192:735-759.
• Bhattacharya, M.C., (1986). Static and dynamic deflections of plates of arbitrary geometry by a new finite difference approach. Journal of Sound and Vibration 107:507-521.
• Chung, C. K., (1981). A generalized finite difference method for heat transfer problems of irregular geometries. Numerical Heat Transfer 4: 345-357.
• Collatz, L., (1966). The Numerical Treatment of Differential Equations. Springer (Berlin).
• Fan. C., Chu. C., Šarler, B., Li. T., (2018). Numerical solutions of waves-current interactions by generalized finite difference method. Engineering Analysis with Boundary Elements. In press: https://doi.org/10.1016/j.enganabound.2018.01.010
  » https://doi.org/10.1016/j.enganabound.2018.01.010
• Fielding, L. M., Villaça, L.M., Garcia, L.F.T., (1997). Energetic finite differences with arbitrary meshes applied to plate-bending problems. Applied Mathematical Modelling 21:691-698.
• Forsythe, G., Wasow., W., (1960). Finite-difference Methods for Partial Differential Equations, John Wiley (New York).
• Gavete, L., Ureña, F., Benito, J.J., Ureña, M., Gavete, M.L., (2018), Solving Elliptical Equations in 3D by Means of an Adaptive Refinement in Generalized Finite Differences. Mathematical Problems in Engineering 2018: Article ID 9678473.
• Ghali, A., Bathe, K. J., (1970). Analysis of plates subjected to in-plane forces using large finite elements. I.A.B.S.E 30:61-72
• Guang-yao, L., Han-bin, Z., (1993). A finite difference method at arbitrary meshes for the bending of plates with variable thickness. Applied Mathematics and Mechanics 14:299-304.
• Huber, M. T., (1923). “Die Theorie der kreuzweise bewehrten Eisenbetonplatten,” Der Bauingenieur, 4, 354-392.
• Jaśkowiec, J., Milewski, S., (2016). Coupling finite element method with meshless finite difference method in thermomechanical problems. Computers and Mathematics with Applications 24:31-155.
• Jensen, P.S. (1972). Finite Difference Techniques for variable grids. Computers & Structures 2:17-29.
• Kao, R., Perrone., N., (1972). Large deflections of flat arbitrary membranes. Computers & Structures 2:535-546.
• Kroon, D., (2011). 2D Line Curvature and Normals. Mathworks File Exchange (Massachusetts).
• Liszka, T., Orkisz, J., (1980). The finite difference method at arbitrary irregular grids and its application in applied mechanics. Computers & Structures 2:83-95.
• Li, P., Fan., C., (2017). Generalized finite difference method for two-dimensional shallow water equations. Engineering Analysis with Boundary Elements 80:58-71.
• Liu, G. R., (2010). Meshfree methods: moving beyond the finite element method. CRC Press (New York), 2nd ed.
• Marcus, H., (1932). Die Theorie elastischer Gewebe und ihre Anwendung auf die Berechnung biegsamer Platten, Springer-Verlag and Heidelberg GmbH and Co. KG (Berlin).
• Pavlin, V., Perrone, N., (1979). Finite difference energy techniques for arbitrary meshes applied to linear plate problems. International Journal for Numerical Methods in Engineering 14:647-664.
• Perrone, N., Kao, R., (1975). A general finite difference method for arbitrary meshes. Computers & Structures 5:45-58.
• Pister, K. S., Dong, S. B., (1959). Elastic Bending of Layered Plates. Journal of Engineering Mechanics Division 85: 1-10.
• Prieto, U. F., Benito, J., Gavete, L., (2011). Application of the generalized finite difference method to solve the advection-diffusion equation. Journal of Computational and Applied Mathematics 235(7):1849-1855.
• Pulino Filho, A. R., (1989). Diferenças Finitas Para Malhas Arbitrárias, Thesis (in Portuguese), University of Campinas, Brazil.
• Srinivasa, N.A., (2006). Adaptive Mesh Refinement for a Finite Difference Scheme Using a Quadtree Decomposition Approach, Texas A&M University, United States.
• Szilard. R., (2004). Theories and Applications of Plate Analysis. John Wiley & Sons, INC (New Jersey).
• Timoshenko, S., Woinowsky-Krieger, S. (1959). Theory of plates and shells, McGraw-Hill (New York).
• Tseng, A. A., Gu, X. S. (1989). A finite difference scheme with arbitrary mesh systems for solving high-order partial differential equations. Computers & Structures 31:319-328.
• Ugural, A. C., (1981). Stresses in Plates and Shells, McGraw-Hill (Boston).
• Ureña, F., Salete. E., Benito. J. J., Gavete, L., (2011). Solving third and fourth order partial differential equations using GFDM. Application to solve problems of plates. International Journal of Computer Mathematics 89:366-376.
• Wei, J., Wang, S., Hou, Q., Dang, J., (2015). Generalized Finite Difference Time Domain Method and Its Application to Acoustics. Mathematical Problems in Engineering 2015: Article ID 640305.
• Zhang, Q., Li, D., Zhang, C., Xu, D., (2016). Multistep finite difference schemes for the variable coefficient delay parabolic equations. Journal of Difference Equations and Applications 22: 745:765.

• Available online: November 05, 2018.

Publication Dates

• Publication in this collection
  2019

History

• Received
  17 July 2018
• Reviewed
  25 Oct 2018
• Accepted
  01 Nov 2018