Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 9142193, 13 pages
https://doi.org/10.1155/2020/9142193

Research Article

Development of High-Order Infinite Element Method for Bending Analysis of Mindlin–Reissner Plates

D. S. Liu,¹ Y. W. Chen,¹ and C. J. Lu²

¹Department of Mechanical Engineering and Advanced Institute of Manufacturing with High-tech Innovations, National Chung Cheng University, Chiayi, Taiwan
²Department of Intelligent Manufacturing Engineering, Guangdong-Taiwan College of Industrial Science & Technology, Dongguan University of Technology, Dongguan, Guangdong, China

Correspondence should be addressed to Y. W. Chen; ywc0314@gmail.com

Received 11 June 2020; Revised 7 September 2020; Accepted 15 September 2020; Published 1 October 2020

Academic Editor: Nhon Nguyen Thanh

Copyright © 2020 D. S. Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

An approach is presented for solving plate bending problems using a high-order infinite element method (IEM) based on Mindlin–Reissner plate theory. In the proposed approach, the computational domain is partitioned into multiple layers of geometrically similar virtual elements which use only the data of the boundary nodes. Based on the similarity, a reduction process is developed to eliminate virtual elements and overcome the problem that the conventional reduction process cannot be directly applied. Several examples of plate bending problems with complicated geometries are reported to illustrate the applicability of the proposed approach and the results are compared with those obtained using ABAQUS software. Finally, the bending behavior of a rectangular plate with a central crack is analyzed to demonstrate that the stress intensity factor (SIF) obtained using the high-order PIEM converges faster and closer than low-order PIEM to the analytical solution.

1. Introduction

The finite element method (FEM) is the most commonly used numerical approach to accurately predict the static and dynamic behaviors of plate structures. The FEM is flexible in identifying solutions to engineering problems that involve complex geometric shapes, different material compositions, and different load forms; therefore, it has become the most extensively implemented numerical method in commercial analysis and simulation software in the market. Nevertheless, the FEM necessitates first establishing an analysis grid containing elements and nodes, a process that is time consuming and labor intensive. Furthermore, to obtain accurate analysis results, a relatively high number of elements and nodes must be established in the calculation domain, and this process exerts a considerable burden on computer memory and reduces the calculation speed. Moreover, direct application of the FEM to the Mindlin–Reissner theory of beams and plates may experience the same numerical error, known as transverse shear locking, that is frequently encountered in FEM analyses. Nguyen-Xuan et al. [1, 2] introduced edge-based and node-based smoothed stabilized discrete shear gap method (ES-DSG, NS-DSG) in conjunction with the high-order shear deformation theory (HSDT) to investigate statics and free vibration behavior of plates; the numerical examples illustrated that both methods are free of shear locking and the results are extremely efficient and accurate.

In the past years, meshless methods have been used [3, 4] as alternatives to the FEM. Such methods require only spatial nodes, obviating the necessity of establishing an analysis grid. Therefore, the disadvantages of the FEM can be overcome, and the analysis time and labor required for extensive modeling tasks can be reduced. Recently, naturally stabilized nodal integration (NSNI) mesh-free formulation has been extensively developed by Thai et al. [5, 6] and successfully used in many complex plate structures such as laminated composite, sandwich plate, and multilayer functionally graded graphene platelets reinforced composite plates. Numerical results show the current approach is promising and highly accurate. Isogeometric analysis (IGA) is a recently introduced technique in the fields of numerical analysis; IGA was first proposed by Hughes et al. [7] as a novel technology to bridge computer-aided design (CAD) and finite element analysis (FEA). Its essential idea is to adopt the same basis functions that are used in geometric design, such as B-splines and nonuniform rational B-splines (NURBS) for the FEM simulations. The combined concept of IGA allows for improved convergence and smoothness properties of the FEM solutions and faster overall simulations. Thus, IGA has been successfully applied to solve demanding problems as geometrically nonlinear analysis [8], buckling and free vibration analysis problem for laminated composite plates [9], and crack growth analysis in thin-shell structures by isogeometric mesh-free coupling approach with a local adaptive mesh refinement scheme near the crack tip [10, 11].

An alternative numerical method called infinite element method (IEM) is a meshless method based on the FEM. In this method, the special similarity between elements can be used to easily create lots of elements as required, and back substitution can be applied to degenerate an infinite number of elements into a multinode super element. Therefore, the IEM can effectively prevent the problems of considerable memory usage and low computing efficiency and speed. The presented method is equally well suited for the usual regularity closed domain and other types of singularities. Furthermore, it can be easily combined with FEM. Thatcher [12, 13] has combined the concept of the FEM and similar splitting to create many tiny elements near a singularity point to approximate Laplace’s equation near a boundary singularity. Moreover, to resolve the problem of structural cracks, Ying and Han [14, 15] have produced many similar triangular elements near a crack tip and combined them into a single element. The calculation results and theoretical solution regarding the stress intensity factor (SIF) were comparable. To solve two-dimensional (2D) and three-dimensional (3D) crack problems, Go et al. [16, 17] have used the similarity of quadrilateral elements to generate so-called super elements by using iterative methods. Liu et al. [18, 19] have combined the IEM with the FEM to solve static linear problems and have continuously extended equations from 2D to 3D. Liu et al. [20] further derived a high-order IEM equation for analyzing various 2D elastic static problems; they compared their results with those of the traditional low-order IEM and with analytical solutions provided in the literature. Their findings revealed that the results obtained using their method were more accurate than those obtained using the low-order IEM and were in good agreement with the analytical solutions. Furthermore, Liu et al. [21] combined the IEM with Mindlin–Reissner plate theory and a closed mode of the IEM to analyze the effects of the size, position, and shape of a circular hole on the flexural stiffness of a thin plate.

Mindlin–Reissner plate theory can be applied to appropriately reduce 3D problems to 2D problems and can be used to increase computing efficiency and reduce memory usage. Currently, this theory is extensively used by scholars. To increase the accuracy and speed of numerical analyses, several scholars have focused on the development of higher order thin plate elements [22, 23]. High-order IEM and Mindlin–Reissner plate theory are the available methods, respectively. However, the conventional reduction process cannot be directly applied when these two theories are combined. Accordingly, a new reduction process has been developed to eliminate virtual elements in the IEM domain so that the IE range is condensed and transformed to form a super element with the master nodes on the boundary only. To demonstrate the effectiveness of the proposed method, we compared the results with that obtained using ABAQUS software. Finally, the analysis results were compared with those obtained using the traditional low-order IEM.

2. Mindlin–Reissner Plate Theory

Mindlin–Reissner plate theory is an extension of Kirchhoff–Love plate theory, which considers shear deformations through the thickness of a plate. When Mindlin–Reissner theory is applied, the following assumptions are used: (a) the thickness of the plate remains unchanged during deformation; (b) the normal stress through the thickness can be ignored; and (c) the normal line of the thickness is perpendicular to the neutral axis line after deformation.

On the basis of the aforementioned assumptions, a complete 3D solid mechanics problem can be reduced to a 2D problem. Therefore, in-plane displacements can be expressed in equations (1) and (2), and the transverse displacement can be expressed as indicated in equation (3).

where x and y are the in-plane axes located in the midplane of the plate and z is the in-plane axis located along the direction of plate thickness (Figure 1). and are the rotations of the midplane about the y and x axes, respectively; and is the angle caused by transverse shear deformation. Executing a transformation from physical to natural coordinates yields the rotation and transverse displacements as follows:

where represents the n-node plate finite element shape function and represents the natural coordinates. The nine-node-plate finite element stiffness matrix can be derived using Mindlin–Reissner theory and by transforming physical coordinates to natural coordinates. The associated plate stiffness is expressed in equation (5), where and denote the bending stiffness and shear stiffness, respectively. The plate material is considered linear elastic, isotropic, and homogenous. The resultant equation of each element can be expressed in equation (12).

where

where h is the plate thickness, κ is the shear energy correction factor (usually 5/6), and is the Jacobian matrix:

Figure 1: Free-body diagram of a plate element.

and comprise shape functions, as presented in equations (8) and (9), respectively. In addition, and are related to the material properties of the model, as presented in equations (10) and (11), respectively.

3. High-Order PIEM

3.1. Similarity Characteristic

Figure 2 presents the basic concept of the infinite element (IE) model. In this model, the computational domain is partitioned into multiple layers of geometrically similar elements. For element I, the local nodes i are numbered 1, 2, …, and 9. If the global origin O and are considered the center of the similarity and the proportionality ratio, respectively, then element II can be created. The global coordinates of elements I and II are related, as presented in equation (13). According to equations (13) and (7), the determinants of the Jacobian matrices of elements I and II are related, as expressed in equation (13). Similarly, according to equations (13) and (8), the relation between of element I and of element II can be presented in equation (15).

Figure 2: Geometrically similar 2D elements in IE formulation.

Therefore, as shown in equation (16), the bending stiffness matrix of the first and second element layers is related.

To adapt the conventional IEM to Mindlin–Reissner plate problems, the shear stiffness of the first element layer can be partitioned into two submatrices, namely, and :

where

Substituting equation (17) into equation (9) yields the following:

Let

Thus, equation (19) becomes

According to the geometric similarity, the relationship between the first and second element layers in terms of the shear stiffness matrix can be expressed as follows:

Substituting equations (16) and (24) into equation (5) yields the plate stiffness matrix as follows:

3.2. Combined Stiffness of High-Order PIEM

According to equation (25), the nine-node elements I, II, …, and s can be mapped using the same square-shaped master element. Specifically, these elements can be designated as similar elements when the coordinate of an element is similar to that of other elements. The matrices of the first element layer can be expressed as follows:

When equation (26) is substituted into equation (12) and the result is expanded, the equations of the s element layers in the computational domain can be derived as follows:

where

In equations (27)–(29), is the nodal displacement vector associated with the ith node layer and is the corresponding nodal force vector. Combining the equations from the first element layer to the sth element layer and assuming that no internal force is applied to the ith node layer (i.e., ) can yield the following expression:

where . If in the last row of equation (31), then

Substituting equation (32) into the second-to-last row of equation (31) yields

Similarly, substituting equations (32) and (33) into the second-to-last row of equation (31) yields

According to equations (32)–(34), the following iteration formulas can be inferred:

Because is equal to , we can iterate , , …, and by using equation (37). According to equation (39), . Substituting into the first row of equation (31) yields the following fundamental IEM formula:

where is the combined stiffness matrix, which preserves the inherent symmetry characteristic of the global stiffness matrix used in the conventional finite element procedure. Using equations (38) and (39), we can condense all inner layer elements and transform them into a single super element with master nodes at the outer boundary.

Ying [24] proved that converges toward a certain constant quantity as the number of element layers approaches infinity; that is,

where s denotes the number of the defined element layers. However, equation (41) cannot be directly applied to the numerical formulation because the infinity element layers are not countable in a physical sense. Therefore, Liu [23] proposed a convergence method for observing the diagonal trace term . The desired accuracy criterion can be expressed as follows:

When this criterion is satisfied, the iterative program is terminated and the critical number of element layers (scr) is determined as equal to the terminated iterative value (i). scr is the minimum number of element layers required for convergence; this implies that sufficient elements are available to cover the entire domain. The proportionality ratio ξ is another important factor in the convergence study. A higher ξ indicates that a higher number of element layers scr is required. Specifically, given a sufficiently high s value, the stiffness is approximately equal to the combined stiffness .

4. Case Studies

4.1. Circular Plate Subject to a Concentrated Load

Consider, for example, a simply supported circular plate subjected to a concentrated load P of 1 lbf at its centroid (Figure 3(a)). The material and geometric parameters are as follows: Young’s modulus E = 3 × 10⁶ psi; Poisson’s ratio ν = 0.3; plate radius R = 10 in; and thickness h = 0.2 in. The analytical maximum deflection was provided by a previous study [25]:

where

(a)

(b)

Figure 3: Circular plate subjected to a point load at the centroid. (a) Analysis model and (b) high-order PIEM.

The solution procedure of the high-order PIEM entails the assumption that the outer boundary comprises 30 uniformly distributed master nodes; in addition, the proportionality ratio ξ is set to 0.64 (Figure 3(b)). On the basis of the convergence criterion (equation (42)), the number of virtual element layers s required is 19. Table 1 illustrates the convergence process. Given the geometric symmetry and load, only a quarter of the entire strip, under the provided load and boundary conditions, must be considered. For comparison, we determine the maximum deflection using ABAQUS software (S4R, 394 elements). The results obtained using ABAQUS are in good agreement with those obtained using the proposed method, as presented in Table 2.

Table 1: Convergence of the required virtual element layer.

Number of element layers
15
16
17
18
19

Accuracy criterion
2.07 × 10⁻⁵
7.81 × 10⁻⁶
2.96 × 10⁻⁶
1.12 × 10⁻⁶
4.22 × 10⁻⁷

Table 2: Maximum deflection of the circular plate subjected to a point load at the centroid.

Methods
Analytical
ABAQUS
High-order PIEM

Maximum deflection
−0.00230
−0.00231
−0.00231

Relative deviation
—
0.49%
0.49%

4.2. Square Plate Subject to a Concentrated Load

Consider a simply supported square plate subjected to a center unit point load (Figure 4(a)). The material and geometric parameters are listed as follows: Young’s modulus E = 3 × 10⁶ psi; Poisson’s ratio ν = 0.3; dimension a = 80 in; and thickness h = 0.8 in. The analytical solution for this problem was provided by a previous study [25], where the deflection at the plate centroid can be expressed as follows:

(a)

(b)

Figure 4: Square plate subjected to a point load at the centroid. (a) Analysis model and (b) high-order PIEM.

In equation (45), the coefficient α (0.0116) is a function of the dimension ratio a = b and D is the flexural rigidity of the plate. The solution procedure of the high-order PIEM involves the assumption that 40 nodes are uniformly distributed and deployed at the boundary; moreover, the proportionality ratio ξ is set to 0.56 (Figure 4(b)). Given the proportionality ratio (0.56), the number of element layers s required to achieve convergence is 33. Because of the geometric symmetry and load, only a quarter of the complete strip, under the given load and boundary conditions, must be considered. For comparison, we also determine the maximum deflection using ABAQUS software (S4R, 400 elements). The results obtained using ABAQUS are in good agreement with the results obtained using the proposed high-order PIEM (Table 3); nevertheless, the results obtained using the high-order PIEM are closer to the analytical solution.

Table 3: Maximum deflection of the square plate subjected to a point load at the centroid.

Methods
Analytical
ABAQUS
High-order PIEM

Maximum deflection
−0.000528
−0.000531
−0.000527

Relative deviation
—
0.56%
0.48%

4.3. Rectangular Plate Subject to Bending Moments

Consider a rectangular plate with the dimensions 100 mm × 50 mm × 0.5 mm (length × width × thickness) (Figure 5(a)). Assume that two of the opposite edges are simply supported and that the other two edges are free such that the applied bending moment (M = 100 N-mm) vanishes along the two simply supported edges. Assume that the plate has a Young’s modulus of 200 GPa and a Poisson’s ratio of 0.3. Figure 5(b) presents the corresponding virtual mesh pattern obtained by the high-order PIEM before the mesh is degenerated to form a single super element. The high-order PIEM solution procedure involves the assumption that the outer boundary comprises 60 uniformly distributed master nodes; furthermore, the proportionality ratio ξ is set to 0.81. On the basis of the convergence criterion (equation (42)), the number of virtual element layers s required is 31. Table 4 presents a comparison of the deflection profile of edge () of the plate (Figure 5(a)) obtained using the proposed high-order PIEM scheme with the profile obtained using ABAQUS. The results obtained from the high-order PIEM are in good agreement with those obtained from ABAQUS (relative deviation < 0.6%) (Table 4).

(a)

(b)

Figure 5: Simply supported rectangular plate subjected to bending moments. (a) Analysis model and (b) high-order PIEM.

Table 4: Comparisons between proposed approach and ABAQUS of maximum deflection along shown in Figure 5.

Coordinate
Deflection (mm)
Relative deviation (%)

High-order PIEM
ABAQUS

−50
0.000
0.000
0.00

−40
0.450
0.452
0.51

−30
0.793
0.795
0.28

−20
1.037
1.039
0.17

−10
1.183
1.184
0.12

0
1.231
1.233
0.11

10
1.183
1.184
0.12

20
1.037
1.039
0.17

30
0.793
0.795
0.28

40
0.450
0.452
0.51

50
0.000
0.000
0.00

4.4. Cantilever Plate Subject to Concentrated Loads

Consider a rectangular plate with the dimensions 100 mm × 50 mm × 0.5 mm (length × width × thickness) (Figure 6(a)). Assume that one of the edges is clamped and that the other three edges are free such that the concentrated loads (P = 5 N) are applied at the end. Moreover, consider that the plate has a Young’s modulus of 200 GPa and a Poisson’s ratio of 0.3. Figure 6(b) presents the corresponding mesh pattern obtained from the high-order PIEM before the mesh is degenerated to form a single super element. The solution procedure of the high-order PIEM involves the assumption that the outer boundary comprises 60 uniformly distributed master nodes; moreover, the proportionality ratio ξ is set to 0.81. Given the proportionality ratio ξ (0.81), the number of element layers s required is 31. Table 5 illustrates a comparison of the deflection profile of edge () of the plate (Figure 6(a)) obtained using the proposed high-order PIEM with the profile obtained using ABAQUS. The two profiles are in good agreement (relative deviation < 0.2%) (Table 5).

(a)

(b)

Figure 6: Cantilever rectangular plate subjected to concentrated loads. (a) Analysis model and (b) high-order PIEM.

Table 5: Comparisons between proposed approach and ABAQUS of maximum deflection along shown in Figure 6.

Coordinate
Deflection (mm)
Relative deviation (%)

High-order PIEM
ABAQUS

−50
0.000
0.000
0.00

−40
−0.311
−0.310
0.12

−30
−1.449
−1.451
0.12

−20
−3.378
−3.383
0.14

−10
−6.000
−6.009
0.15

0
−9.212
−9.225
0.14

10
−12.914
−12.932
0.14

20
−17.008
−17.030
0.13

30
−21.399
−21.425
0.12

40
−25.991
−26.020
0.11

50
−30.689
−30.722
0.11

4.5. L-Shaped Plate Subject to a Concentrated Load

Consider an L-shaped plate with the dimensions 36 mm × 36 mm × 18 mm (L₁ × L₂ × W) (Figure 7(a)), and the thickness is 0.5 mm. The material and geometric parameters are listed as follows: Young’s modulus E = 200000 and Poisson’s ratio ν = 0.3. Assume that two of the opposite edges are clamped, and the concentrated loads (P = 1 N) are applied at the point B. Figure 7(b) presents the corresponding mesh pattern obtained from the high-order PIEM before the mesh is degenerated to form a single super element. The solution procedure of the high-order PIEM involves the assumption that the outer boundary comprises 48 uniformly distributed master nodes. Given the proportionality ratio (0.81), the number of element layers s required to achieve convergence is 34. Table 6 illustrates a comparison of the deflection profile of edge () of the L-shaped plate obtained using the proposed high-order PIEM with the profile obtained using ABAQUS. The two profiles are in good agreement (relative deviation < 0.3%).

(a)

(b)

Figure 7: Cantilever rectangular plate subjected to concentrated loads. (a) Analysis model and (b) high-order PIEM.

Table 6: Comparisons between proposed approach and ABAQUS of maximum deflection along shown in Figure 7.

Coordinate
Deflection (mm)
Relative deviation (%)

High-order PIEM
ABAQUS

−9
−0.1812
−0.1808
0.26

−3
−0.1388
−0.1385
0.20

3
−0.0987
−0.0985
0.15

9
−0.0622
−0.0622
0.11

15
−0.0313
−0.0313
0.06

21
−0.0088
−0.0089
0.28

27
0.0000
0.0000
0.00

4.6. Multihole Plate Subject to Bending Moments

Consider a multihole plate with the dimensions 144 mm × 36 mm × 1 mm (L × W × thickness) (Figure 8), and the circle holes have a radius R = 0.9 mm. The material and geometric parameters are listed as follows: Young’s modulus E = 200000 and Poisson’s ratio ν = 0.3. Assume that two of the opposite edges are simply supported and that the other two edges are free such that the applied bending moment (M = 100 N-mm) vanishes along the two simply supported edges. Figure 9(a) presents the corresponding mesh pattern of one quarter of the complete strip. The model consists of four subdomains, each of which is the IE model. Given the proportionality ratio (0.81), the number of element layers s required to achieve convergence is 34. The one quarter of the complete strip can be a cell and the complete model by combining four cells (Figure 9(b)). Table 7 illustrates a comparison of the deflection profile of edge () of the multihole plate obtained using the proposed high-order PIEM with the profile obtained using ABAQUS. The two profiles are in good agreement (relative deviation < 0.25%). This example not only demonstrates the feasibility of combining IE subdomains, but also copying.

Figure 8: Cantilever rectangular plate subjected to concentrated loads.

(a)

(b)

Figure 9: Mesh pattern of a multihole plate. (a) Quarter model and (b) complete model.

Table 7: Comparisons between the proposed approach and ABAQUS of maximum deflection along shown in Figure 8.

Coordinate
Deflection (mm)
Relative deviation (%)

High-order PIEM
ABAQUS

−72
0.0000
0.0000
0.00

−54
0.2883
0.2877
0.21

−36
0.4813
0.4805
0.16

−18
0.6029
0.6020
0.16

0
0.6388
0.6379
0.15

18
0.6029
0.6020
0.16

36
0.4813
0.4805
0.16

54
0.2883
0.2877
0.21

72
0.0000
0.0000
0.00

4.7. Cracked Plate Subject to Bending Moments

Finally, consider a central crack on a rectangular plate (Figure 10(a)). As boundary conditions, assume that the two edges parallel to the crack are simply supported and moment M is applied to the edges; the other two edges are free. The properties of the plate are as follows: b/h = 2; b/a = 2; c/b = 2; M = 1; Young’s modulus E = 210000; and Poisson’s ratio ν = 0.3. The high- and low-order PIEM algorithms can be used to investigate the stress intensity factor (SIF), which can be evaluated through crack surface displacement extrapolation and can be expressed as follows [26]:

where is the rotations of the midplane about the y axis and r is the distance of the nodal point from the crack tip, as shown in Figure 11, where the stresses near the tip of the crack are modeled using a polar coordinate framework (r, θ) mounted at the crack tip.

(a)

(b)

(c)

Figure 10: Central crack in the rectangular plate. (a) Analysis, (b) high-order, and (c) low-order models.

Figure 11: Modeling of stresses near tip of crack in the elastic body.

The virtual mesh patterns obtained from the high- and low-order PIEM models are displayed in Figures 10(b) and 10(c), respectively. Because of the geometric symmetry and load, only one-half of the complete model must be considered. The outer boundary comprises 101 distributed master nodes. This example uses an open loop model, that is, no elements are generated between the first and last master node, thereby creating a crack. Figure 12 presents the convergence of the high- and low-order PIEM solutions in terms of the required virtual element layers and SIF. The derived results indicate that for more accurate results, the number of element layers and proportionality ratio should be set to 100 and 0.87, respectively, with respect to the crack tip. Table 8 presents the results obtained using the various methods, indicating that the results are in good agreement with those obtained using the proposed method; nevertheless, the results obtained using the high-order PIEM are closer to the analytical solution.

Figure 12: Convergence of high- and low-order PIEM solutions for SIF for a central crack on a plate.

Table 8: SIF for a central crack on a plate subjected to bending moments.

Methods
Ref. [26]
High-order PIEM
Low-order PIEM

SIF
0.9094
0.8895
0.8825

Relative deviation
—
2.19%
2.96%

5. Conclusions

We present a high-order PIEM based on Mindlin–Reissner plate theory. A new reduction process has been developed to eliminate virtual elements in the IEM domain so that the IE range is condensed and transformed to form a super element with the master nodes on the boundary only and overcome the problem that the conventional reduction process cannot be directly applied. Several numerical examples with complex geometries including L-shaped and multihole plates subject to bending moments have been studied; the numerical results are compared with the results obtained using ABAQUS software, and the comparison proves the effectiveness of the proposed scheme. Finally, the bending behavior of a rectangular plate with a central crack is analyzed to demonstrate that the stress intensity factor (SIF) obtained using the high-order PIEM are converge faster and closer to the analytical solution. The numerical results demonstrate that the high-order PIEM provides an accurate and computationally efficient method for analyzing the plate bending problems.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was partially supported by the Advanced Institute of Manufacturing with High-Tech Innovations (AIM-HI) from the Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education (MOE) in Taiwan. This research was also supported by R.O.C. MOST Foundation, Contract nos. MOST109-2634-F-194-004 and MOST109-2634-F-194-001.

References

H. Nguyen-Xuan, G. R. Liu, C. Thai-Hoang, and T. Nguyen-Thoi, “An edge-based smoothed finite element method (ES-FEM) with stabilized discrete shear gap technique for analysis of Reissner-Mindlin plates,” Computer Methods in Applied Mechanics and Engineering, vol. 199, no. 9–12, pp. 471–489, 2010.

C. H. Thai, L. V. Tran, D. T. Tran, T. Nguyen-Thoi, and H. Nguyen-Xuan, “Analysis of laminated composite plates using higher-order shear deformation plate theory and node-based smoothed discrete shear gap method,” Applied Mathematical Modelling, vol. 36, no. 11, pp. 5657–5677, 2012.

T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.

M. Labibzadeh, “Voronoi based discrete least squares meshless method for assessment of stress field in elastic cracked domains,” Journal of Mechanics, vol. 32, no. 3, pp. 267–276, 2016.

C. H. Thai, A. J. M. Ferreira, and H. Nguyen-Xuan, “Naturally stabilized nodal integration meshfree formulations for analysis of laminated composite and sandwich plates,” Composite Structures, vol. 178, pp. 260–276, 2017.

C. H. Thai and P. Phung-Van, “A meshfree approach using naturally atabilized nodal integration for multilayer FG GPLRC complicated plate structures,” Engineering Analysis with Boundary Elements, vol. 117, pp. 346–358, 2020.

T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement,” Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 39–41, pp. 4135–4195, 2005.

T. T. Yu, S. Yin, T. Q. Bui, and S. Hirose, “A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates,” Finite Elements in Analysis and Design, vol. 96, pp. 1–10, 2015.

T. Yu, S. Yin, T. Q. Bui, S. Xia, S. Tanaka, and S. Hirose, “NURBS-based isogeometric analysis of buckling and free vibration problems for laminated composites plates with complicated cutouts using a new simple FSDT theory and level set method,” Thin-Walled Structures, vol. 101, pp. 141–156, 2016.

N. Nguyen‐Thanh, W. Li, and K. Zhou, “Static and free-vibration analyses of cracks in thin-shell structures based on an isogeometric-meshfree coupling approach,” Computational Mechanics, vol. 62, pp. 1287–1309, 2018.

W. Li, N. Nguyen-Thanh, J. Huang, and K. Zhou, “Adaptive analysis of crack propagation in thin-shell structures via an isogeometric-meshfree moving least-squares approach,” Computer Methods in Applied Mechanics and Engineering, vol. 358, 2020.

R. W. Thatcher, “Singularities in the solution of Laplace’s equation in two dimensions,” IMA Journal of Applied Mathematics, vol. 16, no. 3, pp. 303–319, 1975.

R. W. Thatcher, “On the finite element method for unbounded regions,” SIAM Journal on Numerical Analysis, vol. 15, no. 3, pp. 466–477, 1978.

L. A. Ying, “The infinite similar element method for calculating stress intensity factors,” Scientia Sinica, vol. 21, pp. 19–43, 1978.

H. D. Han and L. A. Ying, “An iterative method in the finite element,” Mathematica Numerica Sinica, vol. 1, pp. 91–99, 1979.

C. G. Go and Y. S. Lin, “Infinitely small element for the problem of stress singularity,” Computers & Structures, vol. 37, pp. 547–551, 1991.

C. G. Go and C. C. Guang, “On the use of an infinitely small element for the three-dimensional problem of stress singularity,” Computers & Structures, vol. 45, no. 1, pp. 25–30, 1992.

D. S. Liu and D. Y. Chiou, “A coupled IEM/FEM approach for solving elastic problems with multiple cracks,” International Journal of Solids and Structures, vol. 40, no. 8, pp. 1973–1993, 2003.

D. S. Liu, D. Y. Chiou, and C. H. Lin, “3D IEM formulation with an IEM/FEM coupling scheme for solving elastostatic problems,” Advances in Engineering Software, vol. 34, no. 6, pp. 309–320, 2003.

D.-S. Liu, K.-L. Cheng, and Z.-W. Zhuang, “Development of high-order infinite element method for stress analysis of elastic bodies with singularities,” Journal of Solid Mechanics and Materials Engineering, vol. 4, no. 8, pp. 1131–1146, 2010.

D. S. Liu, C. Y. Tu, and C. L. Chung, “Coupled PIEM/FEM algorithm based on Mindlin-Reissner plate theory for bending analysis of plates with through-thickness hole,” CMES: Computer Modeling in Engineering & Sciences, vol. 92, pp. 573–594, 2013.

H. Nguyen-Xuan, T. Rabczuk, S. Bordas, and J. F. Debongnie, “A smoothed finite element method for plate Analysis,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 13-16, pp. 1184–1203, 2008.

G. R. Liu, X. Y. Cui, and G. Y. Li, “Analysis of mindlin-reissner plates using cell-based smoothed radial point interpolation method,” International Journal of Applied Mechanics, vol. 2, no. 3, pp. 653–680, 2010.

L. A. Ying, Infinite Element Method, Peking University Press, Beijing, China, 1995.

S. P. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, McGraw-Hill, New York, NY, USA, Second edition, 1959.

T. Dirgantara and M. H. Aliabadi, “Stress intensity factors for cracks in thin plates,” Engineering Fracture Mechanics, vol. 69, no. 13, pp. 1465–1486, 2002.