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

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.


Introduction
e finite element method (FEM) is the most commonly used numerical approach to accurately predict the static and dynamic behaviors of plate structures. e 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.
erefore, 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 ai 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. e combined concept of IGA allows for improved convergence and smoothness properties of the FEM solutions and faster overall simulations.
us, 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. erefore, the IEM can effectively prevent the problems of considerable memory usage and low computing efficiency and speed. e 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. atcher [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. e 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 loworder IEM and with analytical solutions provided in the literature. eir 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.

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. erefore, 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). θ x and θ y are the rotations of the midplane about the y and x axes, respectively; and c 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 H i represents the n-node plate finite element shape function and (ζ, η) represents the natural coordinates. e nine-node-plate finite element stiffness matrix can be derived using Mindlin-Reissner theory and by transforming physical coordinates to natural coordinates. e associated plate stiffness is expressed in equation (5), where [K B ] and [K S ] denote the bending stiffness and shear stiffness, respectively. e plate material is considered linear elastic, isotropic, and homogenous. e 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 [J] is the Jacobian matrix: [B B ] and [B S ] comprise shape functions, as presented in equations (8) and (9) (10) and (11), respectively.   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.

Similarity Characteristic.
e 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 [B B ] of element I and [B B ] of element II can be presented in equation (15).
erefore, as shown in equation (16), the bending stiffness matrix [K B ] 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 where Substituting equation (17) into equation (9) yields the following: Let us, 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:

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. e 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), δ i is the nodal displacement vector associated with the ith node layer and f i 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., f s � 0) can yield the following expression: Mathematical Problems in Engineering where Q i � K a(i+1) + K bi . If M s � K bs 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 M s is equal to K bs , we can iterate M s−1 , M s−2 , . . ., and M 1 by using equation (37). According to equation (39), δ 1 � M −1 1 (A 1 δ 0 + C 1 δ m1 ). Substituting δ 1 into the first row of equation (31) yields the following fundamental IEM formula: 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 K z 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. erefore, Liu [23] proposed a convergence method for observing the diagonal trace term K (s) Z . e 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 (s cr ) is determined as equal to the terminated iterative value (i). s cr is the minimum number of element layers required for convergence; this implies that sufficient elements are available to cover the entire domain. e proportionality ratio ξ is another important factor in the convergence study. A higher ξ indicates that a higher number of element layers s cr is required. Specifically, given a sufficiently high s value, the stiffness K (s) Z is approximately equal to the combined stiffness K Z .

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)). e material and geometric parameters are as follows: Young's modulus E � 3 × 10 6 psi; Poisson's ratio ] � 0.3; plate radius R � 10 in; and thickness h � 0.2 in. e analytical maximum deflection was provided by a previous study [25]: where D � Eh 3 e 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). e results obtained using ABAQUS are in good agreement with those obtained using the proposed method, as presented in Table 2.

Square Plate Subject to a Concentrated Load.
Consider a simply supported square plate subjected to a center unit point load (Figure 4(a)). e material and geometric parameters are listed as follows: Young's modulus E � 3 × 10 6 psi; Poisson's ratio ] � 0.3; dimension a � 80 in; and thickness h � 0.8 in. e analytical solution for this problem was provided by a previous study [25], where the deflection at the plate centroid can be expressed as follows: 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. e 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). e 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.

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. e 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 (AB) of the plate (Figure 5(a)) obtained using the proposed highorder PIEM scheme with the profile obtained using ABA-QUS. e results obtained from the high-order PIEM are in good agreement with those obtained from ABAQUS (relative deviation < 0.6%) ( Table 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. e 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 (AB) of the plate (Figure 6(a)) obtained using the proposed highorder PIEM with the profile obtained using ABAQUS. e two profiles are in good agreement (relative deviation < 0.2%) ( Table 5).

L-Shaped Plate Subject to a Concentrated Load.
Consider an L-shaped plate with the dimensions 36 mm × 36 mm × 18 mm (L 1 × L 2 × W) (Figure 7(a)), and the thickness is 0.5 mm. e 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. e 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 (BC) of the L-shaped plate obtained using the proposed high-order PIEM with the profile obtained using ABA-QUS. e two profiles are in good agreement (relative deviation < 0.3%).

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. e 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. e 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. e 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 (AB) of the multihole plate obtained using the proposed high-order PIEM with the profile obtained using ABAQUS. e two profiles are in good agreement (relative deviation < 0.25%).
is example not only demonstrates the feasibility of combining IE subdomains, but also copying.

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. e 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. e 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 θ y 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. e 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. e outer boundary comprises 101 distributed master nodes. is 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. e 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.

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. e numerical results demonstrate that the high-order PIEM provides an accurate and computationally efficient method for analyzing the plate bending problems.

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

12
Mathematical Problems in Engineering