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

Authors: D. S. Liu; Y. W. Chen; C. J. Lu
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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).(1)u=−zθxx,y=−z∂w∂x−γxz,(2)v=−zθyx,y=−z∂w∂y−γyz,(3)w=wx,y,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 γ is the angle caused by transverse shear deformation. Executing a transformation from physical to natural coordinates yields the rotation and transverse displacements as follows:(4)θx=∑i=1nHiζ,ηθxi,θy=∑i=1nHiζ,ηθyi,w=∑i=1nHiζ,ηwi,where Hi 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 KB and KS 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).(5)K=KB+KS,where(6)KB=h312∫−11∫−11BBTDBBBdetJdζdη,KS=κh∫−11∫−11BSTDSBSdetJdζdη,where h is the plate thickness, κ is the shear energy correction factor (usually 5/6), and J is the Jacobian matrix:(7)J=∂x∂ζ∂y∂ζ∂x∂η∂y∂η.Figure 1
Free-body diagram of a plate element.BB and BS comprise shape functions, as presented in equations (8) and (9), respectively. In addition, DB and DS are related to the material properties of the model, as presented in equations (10) and (11), respectively.(8)BB=∂H1∂x00∂H2∂x00⋯∂H9∂x000∂H1∂y00∂H2∂y0⋯0∂H9∂y0∂H1∂y∂H1∂x0∂H2∂y∂H2∂x0⋯∂H9∂y∂H9∂x0,(9)BS=−H10∂H1∂x−H20∂H2∂x⋯−H90∂H9∂x0−H1∂H1∂y0−H2∂H2∂y⋯0−H9∂H9∂y,(10)DB=E1−ν21ν0ν10001−ν2,(11)DS=E21−ν1001,(12)Kθxθyw=MxMyfz.


## 3. High-Order PIEM

### 3.1. Similarity Characteristic

Figure2 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 BB of element I and BB of element II can be presented in equation (15).(13)xiII,yiII=ξxiI,ξyiI,(14)detJII=ξ2detJI,(15)BBII=1ξBBI,Figure 2
Geometrically similar 2D elements in IE formulation.Therefore, as shown in equation (16), the bending stiffness matrix KB of the first and second element layers is related.(16)KBII=h312∫−11∫−11BBIITDBIIBBIIdetJIIdζdη=h312∫−11∫−111ξBBITDBI1ξBBIξ2detJIdζdη=h312∫−11∫−11BBITDBIBBIdetJIdζdη=KBI.To adapt the conventional IEM to Mindlin–Reissner plate problems, the shear stiffness of the first element layerBS can be partitioned into two submatrices, namely, BS∗ and BS∗∗:(17)BS=BS∗+BS∗∗,where(18)Bs∗=00∂H1∂x00∂H2∂x⋯00∂H9∂x00∂H1∂y00∂H2∂y⋯00∂H9∂y,Bs∗=−H100−H200⋯−H9000−H100−H20⋯0−H90.Substituting equation (17) into equation (9) yields the following:(19)KS=κh∫−11∫−11BS∗+BS∗∗TDSBS∗+BS∗∗detJdζdη.Let(20)KS∗=κh∫−11∫−11BS∗TDSBS∗detJdζdη,(21)KS∗∗=κh∫−11∫−11BS∗TDSBS∗∗+BS∗∗TDSBS∗detJdζdη,(22)KS∗∗∗=κh∫−11∫−11BS∗∗TDSBS∗∗detJdζdη.Thus, equation (19) becomes(23)KS=KS∗+KS∗∗+KS∗∗∗.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:(24)KSII=KS∗II+KS∗∗II+KS∗∗∗II=KS∗I+ξKS∗∗I+ξ2KS∗∗∗I.Substituting equations (16) and (24) into equation (5) yields the plate stiffness matrix as follows:(25)KI=KBI+KS∗I+KS∗∗I+KS∗∗∗IKII=KBI+KS∗I+ξKS∗∗I+ξ2KS∗∗∗I⋮Ks=KBI+KS∗I+ξs−1KS∗∗I+ξ2s−1KS∗∗∗I..


### 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:(26)KBI=Ka−BT−AT−BKm−CT−A−CKb,KS∗I=Ka∗−B∗T−A∗T−B∗Km∗−C∗T−A∗−C∗Kb∗,KS∗∗I=Ka∗∗−B∗∗T−A∗∗T−B∗∗Km∗∗−C∗∗T−A∗∗−C∗∗Kb∗∗,KS∗∗∗I=Ka∗∗∗−B∗∗∗T−A∗∗∗T−B∗∗∗Km∗∗∗−C∗∗∗T−A∗∗∗−C∗∗∗Kb∗∗∗.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:(27)Ka1−B1T−A1T−B1Km1−C1T−A1−C1Kb1δ0δm1δ2=f00f2,⋮(28)Kai−BiT−AiT−BiKmi−CiT−Ai−CiKbiδi−1δmiδi=−fi−10fi,⋮(29)Kas−BsT−AsT−BsKms−CsT−As−CsKbsδs−1δmsδs=−fs−10fs,where(30)Kai=Ka+Ka∗+ξi−1Ka∗∗+ξ2i−1Ka∗∗∗,Kbi=Kb+Kb∗+ξi−1Kb∗∗+ξ2i−1Kb∗∗∗,Kmi=Km+Km∗+ξi−1Km∗∗+ξ2i−1Km∗∗∗,Ai=A+A∗+ξi−1A∗∗+ξ2i−1A∗∗∗,Bi=B+B∗+ξi−1B∗∗+ξ2i−1B∗∗∗,Ci=C+C∗+ξi−1C∗∗+ξ2i−1C∗∗∗.In equations (27)–(29), δi is the nodal displacement vector associated with the ith node layer and fi 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., fs=0) can yield the following expression:(31)Ka1−B1T−A1T0000000−B1Km1−C1T0000000−A1−C1Q1−B2T−A2T0000000−B2Km2−C2T0000000−A2−C2Q2−B3T−A3T000⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱00000−Bs−1Kms−1−Cs−1T0000000−As−1−Cs−1Qs−1−BsT−AsT0000000−BsKms−CsT0000000−As−CsKbsδ0δm1δ1δm2δ2⋮δms−1δs−1δmsδs=f00000⋮0000,where Qi=Kai+1+Kbi. If Ms=Kbs in the last row of equation (31), then(32)δs=Ms−1Asδs−1+Csδms.Substituting equation (32) into the second-to-last row of equation (31) yields(33)δms=Kms−CsTMs−1Cs−1Bs+CsTMs−1Asδ4=Mms−1Nsδs−1.Similarly, substituting equations (32) and (33) into the second-to-last row of equation (31) yields(34)δs−1=Qs−1−AsTMs−1As−NsTMms−1Ns−1As−1δs−2+Cs−1δms=Ms−1−1As−1δs−2+Cs−1δms−1.According to equations (32)–(34), the following iteration formulas can be inferred:(35)Mmi=Kmi−CiTMi−1Ci,i=ms,ms−1,…,m1,(36)Ni=Bi+CiTMi−1Ai,i=s,s−1,…,1,(37)Mi=Qi−1−AiTMi−1Ai−NiTMmi−1Ni,i=s−1,s−2,…,1,(38)δmi=Mmi−1Niδi−1,i=s,s−1,…,1,(39)δi=Mi−1Aiδi−1+Ciδmi,i=s,s−1,…,1.BecauseMs is equal to Kbs, we can iterate Ms−1, Ms−2, …, and M1 by using equation (37). According to equation (39), δ1=M1−1A1δ0+C1δm1. Substituting δ1 into the first row of equation (31) yields the following fundamental IEM formula:(40)Ka1−A1TM2−1A1−N1TM1−1N1δ0=Kzδ0=f0,where Kz=Ka1−A1TM2−1A1−N1TM1−1N1 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 Kz converges toward a certain constant quantity as the number of element layers approaches infinity; that is,(41)lims⟶∞Kzs=Kz,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 KZs. The desired accuracy criterion can be expressed as follows:(42)ε=KZi+1−KZiKZi+1×100%≤10−6.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 KZs is approximately equal to the combined stiffness KZ.


## 3.1. Similarity Characteristic

Figure2 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 BB of element I and BB of element II can be presented in equation (15).(13)xiII,yiII=ξxiI,ξyiI,(14)detJII=ξ2detJI,(15)BBII=1ξBBI,Figure 2
Geometrically similar 2D elements in IE formulation.Therefore, as shown in equation (16), the bending stiffness matrix KB of the first and second element layers is related.(16)KBII=h312∫−11∫−11BBIITDBIIBBIIdetJIIdζdη=h312∫−11∫−111ξBBITDBI1ξBBIξ2detJIdζdη=h312∫−11∫−11BBITDBIBBIdetJIdζdη=KBI.To adapt the conventional IEM to Mindlin–Reissner plate problems, the shear stiffness of the first element layerBS can be partitioned into two submatrices, namely, BS∗ and BS∗∗:(17)BS=BS∗+BS∗∗,where(18)Bs∗=00∂H1∂x00∂H2∂x⋯00∂H9∂x00∂H1∂y00∂H2∂y⋯00∂H9∂y,Bs∗=−H100−H200⋯−H9000−H100−H20⋯0−H90.Substituting equation (17) into equation (9) yields the following:(19)KS=κh∫−11∫−11BS∗+BS∗∗TDSBS∗+BS∗∗detJdζdη.Let(20)KS∗=κh∫−11∫−11BS∗TDSBS∗detJdζdη,(21)KS∗∗=κh∫−11∫−11BS∗TDSBS∗∗+BS∗∗TDSBS∗detJdζdη,(22)KS∗∗∗=κh∫−11∫−11BS∗∗TDSBS∗∗detJdζdη.Thus, equation (19) becomes(23)KS=KS∗+KS∗∗+KS∗∗∗.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:(24)KSII=KS∗II+KS∗∗II+KS∗∗∗II=KS∗I+ξKS∗∗I+ξ2KS∗∗∗I.Substituting equations (16) and (24) into equation (5) yields the plate stiffness matrix as follows:(25)KI=KBI+KS∗I+KS∗∗I+KS∗∗∗IKII=KBI+KS∗I+ξKS∗∗I+ξ2KS∗∗∗I⋮Ks=KBI+KS∗I+ξs−1KS∗∗I+ξ2s−1KS∗∗∗I..


## 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:(26)KBI=Ka−BT−AT−BKm−CT−A−CKb,KS∗I=Ka∗−B∗T−A∗T−B∗Km∗−C∗T−A∗−C∗Kb∗,KS∗∗I=Ka∗∗−B∗∗T−A∗∗T−B∗∗Km∗∗−C∗∗T−A∗∗−C∗∗Kb∗∗,KS∗∗∗I=Ka∗∗∗−B∗∗∗T−A∗∗∗T−B∗∗∗Km∗∗∗−C∗∗∗T−A∗∗∗−C∗∗∗Kb∗∗∗.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:(27)Ka1−B1T−A1T−B1Km1−C1T−A1−C1Kb1δ0δm1δ2=f00f2,⋮(28)Kai−BiT−AiT−BiKmi−CiT−Ai−CiKbiδi−1δmiδi=−fi−10fi,⋮(29)Kas−BsT−AsT−BsKms−CsT−As−CsKbsδs−1δmsδs=−fs−10fs,where(30)Kai=Ka+Ka∗+ξi−1Ka∗∗+ξ2i−1Ka∗∗∗,Kbi=Kb+Kb∗+ξi−1Kb∗∗+ξ2i−1Kb∗∗∗,Kmi=Km+Km∗+ξi−1Km∗∗+ξ2i−1Km∗∗∗,Ai=A+A∗+ξi−1A∗∗+ξ2i−1A∗∗∗,Bi=B+B∗+ξi−1B∗∗+ξ2i−1B∗∗∗,Ci=C+C∗+ξi−1C∗∗+ξ2i−1C∗∗∗.In equations (27)–(29), δi is the nodal displacement vector associated with the ith node layer and fi 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., fs=0) can yield the following expression:(31)Ka1−B1T−A1T0000000−B1Km1−C1T0000000−A1−C1Q1−B2T−A2T0000000−B2Km2−C2T0000000−A2−C2Q2−B3T−A3T000⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱00000−Bs−1Kms−1−Cs−1T0000000−As−1−Cs−1Qs−1−BsT−AsT0000000−BsKms−CsT0000000−As−CsKbsδ0δm1δ1δm2δ2⋮δms−1δs−1δmsδs=f00000⋮0000,where Qi=Kai+1+Kbi. If Ms=Kbs in the last row of equation (31), then(32)δs=Ms−1Asδs−1+Csδms.Substituting equation (32) into the second-to-last row of equation (31) yields(33)δms=Kms−CsTMs−1Cs−1Bs+CsTMs−1Asδ4=Mms−1Nsδs−1.Similarly, substituting equations (32) and (33) into the second-to-last row of equation (31) yields(34)δs−1=Qs−1−AsTMs−1As−NsTMms−1Ns−1As−1δs−2+Cs−1δms=Ms−1−1As−1δs−2+Cs−1δms−1.According to equations (32)–(34), the following iteration formulas can be inferred:(35)Mmi=Kmi−CiTMi−1Ci,i=ms,ms−1,…,m1,(36)Ni=Bi+CiTMi−1Ai,i=s,s−1,…,1,(37)Mi=Qi−1−AiTMi−1Ai−NiTMmi−1Ni,i=s−1,s−2,…,1,(38)δmi=Mmi−1Niδi−1,i=s,s−1,…,1,(39)δi=Mi−1Aiδi−1+Ciδmi,i=s,s−1,…,1.BecauseMs is equal to Kbs, we can iterate Ms−1, Ms−2, …, and M1 by using equation (37). According to equation (39), δ1=M1−1A1δ0+C1δm1. Substituting δ1 into the first row of equation (31) yields the following fundamental IEM formula:(40)Ka1−A1TM2−1A1−N1TM1−1N1δ0=Kzδ0=f0,where Kz=Ka1−A1TM2−1A1−N1TM1−1N1 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 Kz converges toward a certain constant quantity as the number of element layers approaches infinity; that is,(41)lims⟶∞Kzs=Kz,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 KZs. The desired accuracy criterion can be expressed as follows:(42)ε=KZi+1−KZiKZi+1×100%≤10−6.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 KZs is approximately equal to the combined stiffness KZ.


## 4. Case Studies

### 4.1. Circular Plate Subject to a Concentrated Load

Consider, for example, a simply supported circular plate subjected to a concentrated loadP of 1 lbf at its centroid (Figure 3(a)). The material and geometric parameters are as follows: Young’s modulus E = 3 × 106 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]:(43)wmax=3+νPR216π1+νD,where(44)D=Eh3121−ν2.Figure 3
Circular plate subjected to a point load at the centroid. (a) Analysis model and (b) high-order PIEM.
(a)(b)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 layers1516171819Accuracy criterionε2.07 × 10−57.81 × 10−62.96 × 10−61.12 × 10−64.22 × 10−7Table 2
Maximum deflection of the circular plate subjected to a point load at the centroid.
MethodsAnalyticalABAQUSHigh-order PIEMMaximum deflectionwmax−0.00230−0.00231−0.00231Relative 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 (Figure4(a)). The material and geometric parameters are listed as follows: Young’s modulus E = 3 × 106 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:(45)wmax=αPa2D.Figure 4
Square plate subjected to a point load at the centroid. (a) Analysis model and (b) high-order PIEM.
(a)(b)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.
MethodsAnalyticalABAQUSHigh-order PIEMMaximum deflectionwmax−0.000528−0.000531−0.000527Relative 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) (Figure5(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 (AB¯) 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).Figure 5
Simply supported rectangular plate subjected to bending moments. (a) Analysis model and (b) high-order PIEM.
(a)(b)Table 4
Comparisons between proposed approach and ABAQUS of maximum deflection alongAB¯ shown in Figure 5.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−500.0000.0000.00−400.4500.4520.51−300.7930.7950.28−201.0371.0390.17−101.1831.1840.1201.2311.2330.11101.1831.1840.12201.0371.0390.17300.7930.7950.28400.4500.4520.51500.0000.0000.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) (Figure6(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 (AB¯) 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).Figure 6
Cantilever rectangular plate subjected to concentrated loads. (a) Analysis model and (b) high-order PIEM.
(a)(b)Table 5
Comparisons between proposed approach and ABAQUS of maximum deflection alongAB¯ shown in Figure 6.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−500.0000.0000.00−40−0.311−0.3100.12−30−1.449−1.4510.12−20−3.378−3.3830.14−10−6.000−6.0090.150−9.212−9.2250.1410−12.914−12.9320.1420−17.008−17.0300.1330−21.399−21.4250.1240−25.991−26.0200.1150−30.689−30.7220.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 (L1 × L2 × 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 (BC¯) 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%).Figure 7
Cantilever rectangular plate subjected to concentrated loads. (a) Analysis model and (b) high-order PIEM.
(a)(b)Table 6
Comparisons between proposed approach and ABAQUS of maximum deflection alongBC¯ shown in Figure 7.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−9−0.1812−0.18080.26−3−0.1388−0.13850.203−0.0987−0.09850.159−0.0622−0.06220.1115−0.0313−0.03130.0621−0.0088−0.00890.28270.00000.00000.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 (AB¯) 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.Figure 9
Mesh pattern of a multihole plate. (a) Quarter model and (b) complete model.
(a)(b)Table 7
Comparisons between the proposed approach and ABAQUS of maximum deflection alongAB¯ shown in Figure 8.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−720.00000.00000.00−540.28830.28770.21−360.48130.48050.16−180.60290.60200.1600.63880.63790.15180.60290.60200.16360.48130.48050.16540.28830.28770.21720.00000.00000.00


### 4.7. Cracked Plate Subject to Bending Moments

Finally, consider a central crack on a rectangular plate (Figure10(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]:(46)SIF=Eh348limr⟶0π2rθyrπ−θyr−π,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.Figure 10
Central crack in the rectangular plate. (a) Analysis, (b) high-order, and (c) low-order models.
(a)(b)(c)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 Figures10(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.
MethodsRef. [26]High-order PIEMLow-order PIEMSIF0.90940.88950.8825Relative deviation—2.19%2.96%


## 4.1. Circular Plate Subject to a Concentrated Load

Consider, for example, a simply supported circular plate subjected to a concentrated loadP of 1 lbf at its centroid (Figure 3(a)). The material and geometric parameters are as follows: Young’s modulus E = 3 × 106 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]:(43)wmax=3+νPR216π1+νD,where(44)D=Eh3121−ν2.Figure 3
Circular plate subjected to a point load at the centroid. (a) Analysis model and (b) high-order PIEM.
(a)(b)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 layers1516171819Accuracy criterionε2.07 × 10−57.81 × 10−62.96 × 10−61.12 × 10−64.22 × 10−7Table 2
Maximum deflection of the circular plate subjected to a point load at the centroid.
MethodsAnalyticalABAQUSHigh-order PIEMMaximum deflectionwmax−0.00230−0.00231−0.00231Relative 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 (Figure4(a)). The material and geometric parameters are listed as follows: Young’s modulus E = 3 × 106 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:(45)wmax=αPa2D.Figure 4
Square plate subjected to a point load at the centroid. (a) Analysis model and (b) high-order PIEM.
(a)(b)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.
MethodsAnalyticalABAQUSHigh-order PIEMMaximum deflectionwmax−0.000528−0.000531−0.000527Relative 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) (Figure5(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 (AB¯) 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).Figure 5
Simply supported rectangular plate subjected to bending moments. (a) Analysis model and (b) high-order PIEM.
(a)(b)Table 4
Comparisons between proposed approach and ABAQUS of maximum deflection alongAB¯ shown in Figure 5.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−500.0000.0000.00−400.4500.4520.51−300.7930.7950.28−201.0371.0390.17−101.1831.1840.1201.2311.2330.11101.1831.1840.12201.0371.0390.17300.7930.7950.28400.4500.4520.51500.0000.0000.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) (Figure6(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 (AB¯) 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).Figure 6
Cantilever rectangular plate subjected to concentrated loads. (a) Analysis model and (b) high-order PIEM.
(a)(b)Table 5
Comparisons between proposed approach and ABAQUS of maximum deflection alongAB¯ shown in Figure 6.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−500.0000.0000.00−40−0.311−0.3100.12−30−1.449−1.4510.12−20−3.378−3.3830.14−10−6.000−6.0090.150−9.212−9.2250.1410−12.914−12.9320.1420−17.008−17.0300.1330−21.399−21.4250.1240−25.991−26.0200.1150−30.689−30.7220.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 (L1 × L2 × 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 (BC¯) 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%).Figure 7
Cantilever rectangular plate subjected to concentrated loads. (a) Analysis model and (b) high-order PIEM.
(a)(b)Table 6
Comparisons between proposed approach and ABAQUS of maximum deflection alongBC¯ shown in Figure 7.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−9−0.1812−0.18080.26−3−0.1388−0.13850.203−0.0987−0.09850.159−0.0622−0.06220.1115−0.0313−0.03130.0621−0.0088−0.00890.28270.00000.00000.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 (AB¯) 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.Figure 9
Mesh pattern of a multihole plate. (a) Quarter model and (b) complete model.
(a)(b)Table 7
Comparisons between the proposed approach and ABAQUS of maximum deflection alongAB¯ shown in Figure 8.
CoordinateDeflection (mm)Relative deviation (%)High-order PIEMABAQUS−720.00000.00000.00−540.28830.28770.21−360.48130.48050.16−180.60290.60200.1600.63880.63790.15180.60290.60200.16360.48130.48050.16540.28830.28770.21720.00000.00000.00


## 4.7. Cracked Plate Subject to Bending Moments

Finally, consider a central crack on a rectangular plate (Figure10(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]:(46)SIF=Eh348limr⟶0π2rθyrπ−θyr−π,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.Figure 10
Central crack in the rectangular plate. (a) Analysis, (b) high-order, and (c) low-order models.
(a)(b)(c)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 Figures10(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.
MethodsRef. [26]High-order PIEMLow-order PIEMSIF0.90940.88950.8825Relative 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.



---

*Source: 9142193-2020-10-01.xml*