# Bending of Nonconforming Thin Plates Based on the Mixed-Order Manifold Method with Background Cells for Integration

Authors: Xin Qu; Lijun Su; Zhijun Liu; Xingqian Xu; Fangfang Diao; Wei Li
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

As it is very difficult to construct conforming plate elements and the solutions achieved with conforming elements yield inferior accuracy to those achieved with nonconforming elements on many occasions, nonconforming elements, especially Adini’s element (ACM element), are often recommended for practical usage. However, the convergence, good numerical accuracy, and high computing efficiency of ACM element with irregular physical boundaries cannot be achieved using either the finite element method (FEM) or the numerical manifold method (NMM). The mixed-order NMM with background cells for integration was developed to analyze the bending of nonconforming thin plates with irregular physical boundaries. Regular meshes were selected to improve the convergence performance; background cells were used to improve the integration accuracy without increasing the degrees of freedom, retaining the efficiency as well; the mixed-order local displacement function was taken to improve the interpolation accuracy. With the penalized formulation fitted to the NMM for Kirchhoff’s thin plate bending, a new scheme was proposed to deal with irregular domain boundaries. Based on the present computational framework, comparisons with other studies were performed by taking several typical examples. The results indicated that the solutions achieved with the proposed NMM rapidly converged to the analytical solutions and their accuracy was vastly superior to that achieved with the FEM and the traditional NMM.


---

## Body

## 1. Introduction

Thin plate, an important engineering structure, is widely employed in various industries, and the bending of such thin plate, being a fourth-order partial differential problem, has become a popular research topic [1–6]. The classical thin plate theory based on Kirchhoff’s assumption demands that both the transverse displacement and normal rotation be continuous, i.e., C1 continuous [7, 8]. This requirement makes it difficult to construct conforming elements. For instance, only the polynomials with an order not less than 5 are able to meet the requirements of continuity even in a general two-dimensional triangular element [9–11]. Furthermore, the solutions achieved with conforming elements yield inferior accuracy to those achieved with nonconforming elements on many occasions [12]. Thus, such nonconforming elements are often recommended for practical usage.Many well-known nonconforming elements have appeared in the history of the FEM, such as the Zienkiewicz element [13], Morley element [14], Veubeke-1 element [15], Veubeke-2 element [15], and ACM element [16]. Such nonconforming elements have realized their corresponding functions in some fields. Unfortunately, the convergence of these nonconforming elements significantly depends on the mesh regularity [17, 18]. For example, the convergence of the Zienkiewicz element can be achieved only if the triangular meshes are parallel to three fixed directions [18, 19]. The convergence of the ACM element depends on the special geometric properties of rectangular meshes [14]. Furthermore, the regular meshes cannot be available for the finite element method (FEM) in some cases, such as in the domain with irregular physical boundaries. In order to overcome the mesh-dependence phenomenon of nonconforming elements, many transformed elements, such as the quasi-conforming element [20] and generalized conforming element [21], were put forward. Although the convergence of these transformed elements can be obtained easily [22–24], the traditional method of constructing displacement functions was not fitted to these transformed elements and good numerical performance cannot be achieved with one element in any situation through numerical experience [25]. An alternative way to solve the convergence problem of nonconforming elements was using the numerical manifold method (NMM) developed by Shi (1991) [26–28], which had been verified by the authors [25, 29–32] taking several classic cases. However, the numerical accuracy of results obtained with nonconforming elements was not satisfactory when we dealt with complex irregular physical boundaries by using the traditional NMM [27].ACM element, a rectangular finite element of classC0, is a well-known nonconforming element with 12 degrees of freedom [33]. And it passes the patch test for polynomials with a degree less than or equal to 2, whatever are the dimensions of the rectangles when all rectangles are equal [18]. Compared with conforming rectangular elements, ACM element has a lower order shape function, less number of degrees of freedom, and higher accuracy. However, the convergence, good numerical accuracy, and high computing efficiency of ACM element with irregular physical boundaries cannot be achieved using either the FEM or the traditional NMM.In order to solve the above problems, the mixed-order NMM with background cells for integration was developed to analyze the bending of nonconforming thin plates with irregular physical boundaries. Compared with the FEM, the regular meshes can always be selected to improve the convergence performance since two completely independent cover systems are used in the NMM. Different from the traditionally local refinement technology, background cells can be used to improve the integration accuracy without increasing the degrees of freedom, retaining efficiency as well. In contrast to the traditional NMM, the approximation to the real boundary and mesh partitioning in preprocessing are much simpler in the mixed-order NMM with background cells for integration. Compared with the traditional method of upgrading the order of the local displacement function, the mixed-order upgrading method based on Taylor expansion is more flexible and efficient. Then, the first-order approximation was taken to construct a high-order local displacement function to improve the interpolation accuracy of the irregular domain, while the zero-order approximation was employed to retain efficiency in the regular domain. In this way, the best numerical accuracy was achieved, and the convergence and high computing efficiency were correspondingly obtained. With the penalized variational formulation fitted to the NMM for Kirchhoff’s thin plate problems, a new scheme was proposed to deal with irregular domain boundaries and the corresponding computational procedure was introduced in this paper. Finally, several typical examples were presented to verify the mixed-order NMM with background cells for integration by comparing with the results achieved using the FEM and the traditional NMM.


## 2. The Penalized Variational Formulation

In general, the integration over the entire problem is divided into a summation of integrations over all the elements, and the solution is constructed by an interpolation polynomial. However, the continuity of the normal rotation on element surfaces in a nonconforming thin plate cannot be achieved. Thus, the continuous requirements for the variation of the interpolation function cannot be satisfied in this case, which leads to a linearδ-type singular distribution of some derivatives. Importantly, the legality of using the Green theorem cannot be guaranteed, and the final numerical solutions may not always converge to the analytical solutions.The “patch test” introduced by Bazeley et al. (1965) [13] is considered as a standard approach for testing the convergence of nonconforming elements. If an element can pass the “patch test”, then, its convergence can be achieved [34]. The existing studies [12, 17, 35] and numerical practices have shown that the elements formed by regular meshes, as proposed by the NMM, are able to pass the “patch test”. However, the regular meshes cannot be available for the FEM in some cases, particularly in the domain with irregular physical boundaries, which creates convergence issues for the nonconforming elements. Since two completely independent cover systems are used in the NMM, the regular meshes can always be selected as the mathematical meshes for interpolation to improve the convergence performance. Even though the physical and mathematical meshes are mismatched, the best meshes are still available after the minimum potential energy principle of the solution domain is slightly modified.Based on Kirchhoff’s thin plate theory [36], the equilibrium equation in terms of w (x, y) under the coordinate system shown in Figure 1 can be given as follows:(1)D∇2w=p¯.Here,(2)D=Eh3121−μ2,(3)∇2=∂x2+∂y22,where w is the deflection function, p¯ is the uniformly distributed load of intensity, D is the flexural rigidity of the plate, h is the thickness of the plate, E is the modulus of elasticity, and μ is Poisson’s ratio.Figure 1
Coordinate system and internal forces.Figure 2
Definition of local frame.The local frame on the boundary was constituted byn and s and displayed in Figure 2. The exterior unit normal vector n and the unit tangential vector s are defined as follows:(4)n=nx,nyΤ,s=−ny,nxΤ.The transformations of first-order derivatives can be written as follows:(5)∂∂n=nx∂∂x+ny∂∂y,∂∂s=−ny∂∂x+nx∂∂y.In general, the boundary of Ω can be divided into three disconnected parts; that is, ∂Ω= S1∪S2∪S3 [28]. S1 is clamped, if on S1: (6)w=w¯,∂w∂n=φ¯n.Here,w¯ is a deflection and φ¯n is a normal angle. S2 is simply supported, if(7)w=w¯,Mn=M¯n.Here,M¯n is a bending moment. S3 is free, if on S3: (8)Mn=M¯n,∂Mns∂s+Qn=q¯.Here,Mns is a twisting moment, Qn is a lateral shear force, and q¯ is a linearly transverse load. Additionally, a quantity under the bar ‘–’ indicates that it is known.Thus, the total potential energy based on Kirchhoff’s thin plate theory can be given as follows:(9)∏p=12∬ΩκTDκ−p¯wdxdy−∫S3q¯wdS+∫SMM¯n∂w∂ndS,where κ is the generalized strain, and SM = S2∪S3Here,(10)κ=−∂2w∂x2−∂2w∂y2−2∂2w∂x∂y,D=D1μ0μ1001−μ.When imposing essential boundary conditions, the FEM requires that all nodes must be on the domain boundaries. This creates difficulties when dealing with the complex boundary, because the considerable effort is required to solve the mismatch between the meshes and domain boundaries. Fortunately, by using the Lagrange multiplier method and penalty method, the essential boundary conditions can be weakly imposed in the NMM [25]. In this paper, the essential boundary conditions were imposed by using the penalty method. The penalized formulation fitted to the total potential energy can be given as follows:(11)∏p∗=∏p+∫Sw12kww−w¯2dS+∫S112kM∂w∂n−φn¯2dS,where Sw = S1∪S2 is the essential boundary and kw and kM are the penalty parameters, varying little for the range of 103E–105E.


## 3. Element Analysis

A series of the identical rectangular meshes were employed as mathematical meshes to cover the neutral surface of a thin plate. In the mathematical meshes, one mathematical patch corresponding to one node is the domain composed of rectangular elements connected to it (customarily called a “mathematical star”). Taking Figure3 as an example, “mathematical star” 1 is the large rectangle 7-8-9-2-3-4-5-6 that contains the four small rectangles 7-8-1-6, 8-9-2-1, 1-2-3-4, and 6-1-4-5; “mathematical star” 2 is the large rectangle 8-9-10-11-12-3-4-1 that contains the four small rectangles 8-9-2-1, 9-10-11-2, 2-11-12-3, and 1-2-3-4; “mathematical star” 3 is the large rectangle 1-2-11-12-13-14-15-4 that contains the four small rectangles 1-2-3-4, 2-11-12-3, 3-12-13-14, and 4-3-14-15; and “mathematical star” 4 is the large rectangle 6-1-2-3-14-15-16-5 that contains the four small rectangles 6-1-4-5, 1-2-3-4, 4-3-14-15, and 5-4-15-16.Figure 3
A manifold element and its “physical stars”.By intersecting the neutral surface of the thin plate with all “mathematical stars”, “physical stars” can be achieved. All “physical stars” constitute the physical covers or physical meshes. The overlapping domain of four physical stars forms a manifold element. In Figure3, the overlapping domain of the physical stars 1, 2, 3, and 4 is the rectangle 1-2-3-4, which is a manifold element.By cohering four “physical stars” with their corresponding weighted (shape) functions, the global deflection function of an element can be obtained as follows:(12)w=∑i=14Niδii=1,2,3,4,where Ni is the weighted (shape) function of “physical star” i, and δi is the local displacement function of “physical star” i.


### 3.1. Shape Function

In this paper, the center of a rectangle served as the origin in the local coordinate system, and theξ and η axes were parallel to the x and y axes, respectively, shown in Figure 4. Then, the transformation relations of these two coordinate systems can be given as follows:where (xi, yi) (i = 1, 2, 3, 4) is the global coordinate of node i.Figure 4
Global and local coordinate systems.(13)x=x0+aξ,x0=x1+x32=x2+x42,y0=y1+y32=y2+y42,y=y0+bη,a=x3−x12=x4−x22,b=y3−y12=y4−y22,In the local coordinate system, the coordinate of nodei is (ξi, ηi) (i = 1, 2, 3, 4), and its value is ±1. In this paper, the commonly used shape functions in the local coordinate ξoη were adopted as follows:(14)Ni=1+ξ01+η02+ξ0+η0−ξ2−η28,Nxi=−bηi1+ξ01+η01−η28,i=1,2,3,4,Nyi=aξi1+ξ01+η01−ξ28,where ξ0 = ξξi and η0 = ηηi. And it is easy to prove that the element is a nonconforming element by adopting the abovementioned shape functions [13, 18, 32]. The shape functions with 0.2 × 0.2 meshes are shown in Figure 5.Figure 5
3D view of the shape functions.


### 3.2. Local Displacement Function

#### 3.2.1. Zero-Order Local Displacement Function

In principle, it is available for each “physical star” to be assigned independently a local deflection function,wi (x, y), a local rotation function about the x axis, θxi (x, y), and a local rotation function about the y axis, θyi (x, y). The undetermined constants in these three functions are uniformly characterized as the generalized degrees of freedom. When constants are chosen as the deflection and rotation functions, it renders the so-called zero-order approximation. In most regular domains, good accuracy and efficiency can be retained using the zero-order approximation.For the zero-order approximation, the local displacement functions are three undetermined constants. Thus, every “physical star” has three generalized degrees of freedom. The degrees of freedom array of “physical star”i can be written as follows:(15)δi=wiθxiθyi.Since every manifold element is formed by four overlapping “physical stars”, there are 12 degrees of freedom in a manifold element. Substituting equation (15) into equation (12), the global deflection function can be given as follows:(16)w=∑i=14Niδi=∑i=14Niwi+Nxiθxi+Nyiθyi,where Ni = [Ni, Nxi, Nyi] must meet equation (17) as follows:(17)Nirj=δij,∂Nirj∂x=0,∂Nirj∂y=0,Nxirj=0,∂Nxirj∂x=0,∂Nxirj∂y=δij,Nyirj=0,∂Nyirj∂x=−δij,∂Nyirj∂y=0.The physical relations between the global deflection function and the global rotation function can be expressed as follows:(18)−θy=∂wx,y∂x,θx=∂wx,y∂y,where θx and θy denote the global rotation functions about x and y axes, respectively.Combining equations (16)–(18), the partial derivatives of the deflection function can be achieved as follows:(19)θxi=∂w∂y|i=∂Ni∂ywi+∂Nxi∂yθxi+∂Nyi∂yθyi=θxi,(20)θxi=−∂w∂x|i=−∂Ni∂xwi−∂Nxi∂xθxi−∂Nyi∂xθyi=θyi.Equations (19) and (20) represent the relations between the global rotation function and the local rotation function at node i in the zero-order approximation.


#### 3.2.2. First-Order Local Displacement Function

In the irregular domain, good accuracy can be seldom achieved using the zero-order approximation. Therefore, a high-order local displacement function is encouraged in such domain. When a linear polynomial is selected as the local deflection function and constants are employed as the local rotation functions about thex and y axes, it is called the first-order approximation. Obviously, the first-order approximation is the most economical way to upgrade the order of the local deflection function. Thus, the first-order approximation was selected to construct a high-order local displacement function to improve the interpolation accuracy of the irregular domain.To endow the generalized degrees of freedom with physical meanings, the local deflection function was constructed following a previously reported method [37, 38]. The local deflection function wi (x, y) was regarded as the first-order Taylor expansion of the global deflection function w (x, y) at (xi, yi), which can be written as follows:(21)wix,y=wi+wxix−xi+wyiy−yi,(22)wi=wxi,yi,wxi=∂wx,y∂x|xi,yi,wyi=∂wx,y∂y|xi,yi.Then, equation (2) can be also written as follows:(23)wix,y=wi−θyix−xi+θxiy−yi,where θxi and θyi, induced by the local deflection function wi (x, y) at (xi, yi), are called the additional rotation functions.Therefore, the local displacement function can be written as follows:(24)wix,yθxiθyi=100−x−xiy−yi0100000100wiθxiθyiθyiθxi,where θxi and θyi denote the local rotation functions about the x and y axes, respectively.From equation (24), the generalized degrees of freedom array of “physical star” i in the first-order approximation can be written as follows:(25)δi=wiθxiθyi=wiθxiθyiθyiθxi.Since every manifold element is formed by four overlapping “physical stars”, there are 20 degrees of freedom in a manifold element. Substituting equation (25) into equation (21), the global deflection function of an element can be given as follows:(26)w=∑i=14NiBi=∑i=14Niwi+Nxiθxi+Nyiθyi+Niθxiy−yi−Niθyix−xi.Combining equations (17), (18), and (26), the partial derivatives of deflection function can be obtained:(27)θxi=∂w∂y|i=∂Ni∂ywi+∂Nxi∂yθxi+∂Nyi∂yθyi+Niθxi=θxi+θxi,(28)θyi=−∂w∂x|i=−∂Ni∂xwi−∂Nxi∂xθxi−∂Nyi∂xθyi+Niθyi=θyi+θyi.According to equations (27) and (28), the global rotations (θx) i and (θy) i contained not only the local rotations θxi and θyi but also the additional rotations θxi and θyi. The reason for this was that the partial derivatives of the first-order local deflection function were not zero, and thus, additional terms appeared after differentiating the global deflection function, which corresponded to the additional rotations θxi and θyi. In such a case, the relations between the global and local rotation functions in the zero-order approximation, shown in equations (19) and (20), did not hold anymore.


## 3.1. Shape Function

In this paper, the center of a rectangle served as the origin in the local coordinate system, and theξ and η axes were parallel to the x and y axes, respectively, shown in Figure 4. Then, the transformation relations of these two coordinate systems can be given as follows:where (xi, yi) (i = 1, 2, 3, 4) is the global coordinate of node i.Figure 4
Global and local coordinate systems.(13)x=x0+aξ,x0=x1+x32=x2+x42,y0=y1+y32=y2+y42,y=y0+bη,a=x3−x12=x4−x22,b=y3−y12=y4−y22,In the local coordinate system, the coordinate of nodei is (ξi, ηi) (i = 1, 2, 3, 4), and its value is ±1. In this paper, the commonly used shape functions in the local coordinate ξoη were adopted as follows:(14)Ni=1+ξ01+η02+ξ0+η0−ξ2−η28,Nxi=−bηi1+ξ01+η01−η28,i=1,2,3,4,Nyi=aξi1+ξ01+η01−ξ28,where ξ0 = ξξi and η0 = ηηi. And it is easy to prove that the element is a nonconforming element by adopting the abovementioned shape functions [13, 18, 32]. The shape functions with 0.2 × 0.2 meshes are shown in Figure 5.Figure 5
3D view of the shape functions.


## 3.2. Local Displacement Function

### 3.2.1. Zero-Order Local Displacement Function

In principle, it is available for each “physical star” to be assigned independently a local deflection function,wi (x, y), a local rotation function about the x axis, θxi (x, y), and a local rotation function about the y axis, θyi (x, y). The undetermined constants in these three functions are uniformly characterized as the generalized degrees of freedom. When constants are chosen as the deflection and rotation functions, it renders the so-called zero-order approximation. In most regular domains, good accuracy and efficiency can be retained using the zero-order approximation.For the zero-order approximation, the local displacement functions are three undetermined constants. Thus, every “physical star” has three generalized degrees of freedom. The degrees of freedom array of “physical star”i can be written as follows:(15)δi=wiθxiθyi.Since every manifold element is formed by four overlapping “physical stars”, there are 12 degrees of freedom in a manifold element. Substituting equation (15) into equation (12), the global deflection function can be given as follows:(16)w=∑i=14Niδi=∑i=14Niwi+Nxiθxi+Nyiθyi,where Ni = [Ni, Nxi, Nyi] must meet equation (17) as follows:(17)Nirj=δij,∂Nirj∂x=0,∂Nirj∂y=0,Nxirj=0,∂Nxirj∂x=0,∂Nxirj∂y=δij,Nyirj=0,∂Nyirj∂x=−δij,∂Nyirj∂y=0.The physical relations between the global deflection function and the global rotation function can be expressed as follows:(18)−θy=∂wx,y∂x,θx=∂wx,y∂y,where θx and θy denote the global rotation functions about x and y axes, respectively.Combining equations (16)–(18), the partial derivatives of the deflection function can be achieved as follows:(19)θxi=∂w∂y|i=∂Ni∂ywi+∂Nxi∂yθxi+∂Nyi∂yθyi=θxi,(20)θxi=−∂w∂x|i=−∂Ni∂xwi−∂Nxi∂xθxi−∂Nyi∂xθyi=θyi.Equations (19) and (20) represent the relations between the global rotation function and the local rotation function at node i in the zero-order approximation.


### 3.2.2. First-Order Local Displacement Function

In the irregular domain, good accuracy can be seldom achieved using the zero-order approximation. Therefore, a high-order local displacement function is encouraged in such domain. When a linear polynomial is selected as the local deflection function and constants are employed as the local rotation functions about thex and y axes, it is called the first-order approximation. Obviously, the first-order approximation is the most economical way to upgrade the order of the local deflection function. Thus, the first-order approximation was selected to construct a high-order local displacement function to improve the interpolation accuracy of the irregular domain.To endow the generalized degrees of freedom with physical meanings, the local deflection function was constructed following a previously reported method [37, 38]. The local deflection function wi (x, y) was regarded as the first-order Taylor expansion of the global deflection function w (x, y) at (xi, yi), which can be written as follows:(21)wix,y=wi+wxix−xi+wyiy−yi,(22)wi=wxi,yi,wxi=∂wx,y∂x|xi,yi,wyi=∂wx,y∂y|xi,yi.Then, equation (2) can be also written as follows:(23)wix,y=wi−θyix−xi+θxiy−yi,where θxi and θyi, induced by the local deflection function wi (x, y) at (xi, yi), are called the additional rotation functions.Therefore, the local displacement function can be written as follows:(24)wix,yθxiθyi=100−x−xiy−yi0100000100wiθxiθyiθyiθxi,where θxi and θyi denote the local rotation functions about the x and y axes, respectively.From equation (24), the generalized degrees of freedom array of “physical star” i in the first-order approximation can be written as follows:(25)δi=wiθxiθyi=wiθxiθyiθyiθxi.Since every manifold element is formed by four overlapping “physical stars”, there are 20 degrees of freedom in a manifold element. Substituting equation (25) into equation (21), the global deflection function of an element can be given as follows:(26)w=∑i=14NiBi=∑i=14Niwi+Nxiθxi+Nyiθyi+Niθxiy−yi−Niθyix−xi.Combining equations (17), (18), and (26), the partial derivatives of deflection function can be obtained:(27)θxi=∂w∂y|i=∂Ni∂ywi+∂Nxi∂yθxi+∂Nyi∂yθyi+Niθxi=θxi+θxi,(28)θyi=−∂w∂x|i=−∂Ni∂xwi−∂Nxi∂xθxi−∂Nyi∂xθyi+Niθyi=θyi+θyi.According to equations (27) and (28), the global rotations (θx) i and (θy) i contained not only the local rotations θxi and θyi but also the additional rotations θxi and θyi. The reason for this was that the partial derivatives of the first-order local deflection function were not zero, and thus, additional terms appeared after differentiating the global deflection function, which corresponded to the additional rotations θxi and θyi. In such a case, the relations between the global and local rotation functions in the zero-order approximation, shown in equations (19) and (20), did not hold anymore.


## 3.2.1. Zero-Order Local Displacement Function

In principle, it is available for each “physical star” to be assigned independently a local deflection function,wi (x, y), a local rotation function about the x axis, θxi (x, y), and a local rotation function about the y axis, θyi (x, y). The undetermined constants in these three functions are uniformly characterized as the generalized degrees of freedom. When constants are chosen as the deflection and rotation functions, it renders the so-called zero-order approximation. In most regular domains, good accuracy and efficiency can be retained using the zero-order approximation.For the zero-order approximation, the local displacement functions are three undetermined constants. Thus, every “physical star” has three generalized degrees of freedom. The degrees of freedom array of “physical star”i can be written as follows:(15)δi=wiθxiθyi.Since every manifold element is formed by four overlapping “physical stars”, there are 12 degrees of freedom in a manifold element. Substituting equation (15) into equation (12), the global deflection function can be given as follows:(16)w=∑i=14Niδi=∑i=14Niwi+Nxiθxi+Nyiθyi,where Ni = [Ni, Nxi, Nyi] must meet equation (17) as follows:(17)Nirj=δij,∂Nirj∂x=0,∂Nirj∂y=0,Nxirj=0,∂Nxirj∂x=0,∂Nxirj∂y=δij,Nyirj=0,∂Nyirj∂x=−δij,∂Nyirj∂y=0.The physical relations between the global deflection function and the global rotation function can be expressed as follows:(18)−θy=∂wx,y∂x,θx=∂wx,y∂y,where θx and θy denote the global rotation functions about x and y axes, respectively.Combining equations (16)–(18), the partial derivatives of the deflection function can be achieved as follows:(19)θxi=∂w∂y|i=∂Ni∂ywi+∂Nxi∂yθxi+∂Nyi∂yθyi=θxi,(20)θxi=−∂w∂x|i=−∂Ni∂xwi−∂Nxi∂xθxi−∂Nyi∂xθyi=θyi.Equations (19) and (20) represent the relations between the global rotation function and the local rotation function at node i in the zero-order approximation.


## 3.2.2. First-Order Local Displacement Function

In the irregular domain, good accuracy can be seldom achieved using the zero-order approximation. Therefore, a high-order local displacement function is encouraged in such domain. When a linear polynomial is selected as the local deflection function and constants are employed as the local rotation functions about thex and y axes, it is called the first-order approximation. Obviously, the first-order approximation is the most economical way to upgrade the order of the local deflection function. Thus, the first-order approximation was selected to construct a high-order local displacement function to improve the interpolation accuracy of the irregular domain.To endow the generalized degrees of freedom with physical meanings, the local deflection function was constructed following a previously reported method [37, 38]. The local deflection function wi (x, y) was regarded as the first-order Taylor expansion of the global deflection function w (x, y) at (xi, yi), which can be written as follows:(21)wix,y=wi+wxix−xi+wyiy−yi,(22)wi=wxi,yi,wxi=∂wx,y∂x|xi,yi,wyi=∂wx,y∂y|xi,yi.Then, equation (2) can be also written as follows:(23)wix,y=wi−θyix−xi+θxiy−yi,where θxi and θyi, induced by the local deflection function wi (x, y) at (xi, yi), are called the additional rotation functions.Therefore, the local displacement function can be written as follows:(24)wix,yθxiθyi=100−x−xiy−yi0100000100wiθxiθyiθyiθxi,where θxi and θyi denote the local rotation functions about the x and y axes, respectively.From equation (24), the generalized degrees of freedom array of “physical star” i in the first-order approximation can be written as follows:(25)δi=wiθxiθyi=wiθxiθyiθyiθxi.Since every manifold element is formed by four overlapping “physical stars”, there are 20 degrees of freedom in a manifold element. Substituting equation (25) into equation (21), the global deflection function of an element can be given as follows:(26)w=∑i=14NiBi=∑i=14Niwi+Nxiθxi+Nyiθyi+Niθxiy−yi−Niθyix−xi.Combining equations (17), (18), and (26), the partial derivatives of deflection function can be obtained:(27)θxi=∂w∂y|i=∂Ni∂ywi+∂Nxi∂yθxi+∂Nyi∂yθyi+Niθxi=θxi+θxi,(28)θyi=−∂w∂x|i=−∂Ni∂xwi−∂Nxi∂xθxi−∂Nyi∂xθyi+Niθyi=θyi+θyi.According to equations (27) and (28), the global rotations (θx) i and (θy) i contained not only the local rotations θxi and θyi but also the additional rotations θxi and θyi. The reason for this was that the partial derivatives of the first-order local deflection function were not zero, and thus, additional terms appeared after differentiating the global deflection function, which corresponded to the additional rotations θxi and θyi. In such a case, the relations between the global and local rotation functions in the zero-order approximation, shown in equations (19) and (20), did not hold anymore.


## 4. Discretization Scheme of Penalized Variational Formulation

According to the principle of minimum potential energy, lettingδ∏pw∗=0, the true deflection was obtained. The discretized algebraic equations can be written as follows:(29)Kδ=f,where K is the global stiffness matrix, and f is the equivalent nodal force vector. Matrix K is obtained through assembling the stiffness matrix Ke from all elements and all displacement boundary segments. Vector f is obtained through assembling the vector fe.The nodal momentM can be derived as follows:(30)M=Dκ,where the generalized strain κ can be written as follows:(31)κ=Bδe=B1B2B3B4δ1Τδ2Τδ3Τδ4ΤΤ,B=B1B2B3B4,δe=δ1Τδ2Τδ3Τδ4ΤΤ,Bi=−∂2Ni∂x2∂2Ni∂y22∂2Ni∂x∂y=−1a2∂2Ni∂ξ21b2∂2Ni∂η22ab∂2Ni∂ξ∂η=−1abba∂2Ni∂ξ2ab∂2Ni∂η22∂2Ni∂ξ∂η.From equations (9), (11), and (21), the element stiffness matrix Ke was composed of the system stiffness matrix k, the matrices due to constrained deflection kw, and normal rotation kθ as follows:(32)Ke=k+kw+kθ,where the system stiffness matrix k can be divided into submatrices as follows:(33)k=k11k12k13k14k21k22k23k24k31k32k33k34k41k42k43k44,kij=∫ΩeBiΤDBjdxdy.Thekw and kθ can be expressed as follows:(34)kw=∫SwkwNΤNdS,kθ=∫S1kM∂NΤ∂n∂N∂ndS.The equivalent load vector can be given as follows:(35)fe=Qe+∫Sww¯kwNΤdS+∫Swφn¯kM∂NΤ∂ndS,Qe=Q1TQ2TQ3TQ4TΤ.Suppose that a thin plate is subjected to a bending momentM¯, linearly transverse load q¯, and uniformly distributed load of intensity p¯, and thus, the equivalent nodal forces can be written as follows:(36)Qi=FziMθxiMθyi=∬ΩeNiΤp¯dxdy+∫S3NiΤq¯eds−∫SM∂NiΤ∂nM¯dsi=1,2,3,4.


## 5. Background Cells for Integration

For the domains with irregular physical boundaries, especially with curved boundaries, the traditional meshes are difficult to match, which leads to an inaccurate approximation and deteriorates the accuracy of solutions. Thus, locally refined meshes are often adopted to approximate these irregular physical boundaries.However, the traditionally local refinement technology must generate more degrees of freedom in the FEM, which deteriorates the efficiency. In addition, the problem of hanging nodes needs to be solved, and preprocessing is more complex when the traditionally local refinement technology is adopted in the NMM. Since the NMM creates two completely independent cover systems (the physical cover is used for integration and the mathematical cover is used for interpolation), considerably refined meshes designated as background cells can be used to approximate the domains with irregular physical boundaries, while relatively coarse meshes were employed for interpolation. In this way, the best integration accuracy can be achieved without increasing the degrees of freedom, retaining efficiency. It should be noted that the traditional NMM did not provide a scheme for refined meshes (background cells) generation, and thus, a new mesh generation technology was developed and introduced here.First, the irregular physical boundaries were approximated by a series of straight lines. In this process, we did not need to have a very accurate approximation, which can greatly simplify the preprocessing and shorten the time for mesh partitioning. Although the region composed of the initial approaching boundary and real boundary is neglected in the traditional NMM, leading to inaccurate integration domains, the inaccurate integration domains can be revised by adding background cells. As shown in Figure6, the circular boundary, called the real boundary, represented by the red thick curve, is approximated by a series of blue thick straight lines called the initial approaching boundary.Figure 6
Relations of the initial approaching boundary, the real boundary, and rectangular meshes.Then, a series of rectangular meshes were selected to cover the domains, and the meshes for integration were generated by cutting the rectangular meshes with the initial approaching boundaries. As shown in Figure6, the intersections between the initial approaching boundary and the rectangular meshes were classified into two groups, i.e., A1,A2,…,An (on the real boundary), and B1,B2,…,Bm (not on the real boundary). The nodes B1,B2,…,Bm will be replaced by the intersections between the real boundary and the rectangular meshes, D1,D2,…,Dk, which can improve the accuracy of boundary approximation. The nodes D1,D2,…,Dk were the closest to B1,B2,…,Bm. The nodes C1,C2,…,Cl lied inside the irregular elements, and they were the intersections between the initial approaching boundary and the real boundary. The nodes A1,A2,…,An and C1,C2,…,Cl were sorted anticlockwise, and then they were recorded as E1,E2,…,EN (N = n + l). Every two adjacent nodes, Ei, Ei+1 (1 ≤i≤ N − 1), were connected into a straight line Ei-Ei + 1. The real domain was cut into two parts by straight line Ei-Ei + 1, in which the smaller one was recorded as Ωi. The nodes D1,D2,…,Dk lied in Ωi were figured out and sorted anticlockwise, recorded as {F1, F2,…,FK}. The nodes {F1, F2,…,FK}, Ei, Ei+1 were connected to construct polygon Ei-F1-F2-,…,-FK-Ei + 1. If the ratio of the area of polygon Ei-F1-F2-,…,-FK-Ei+1 to the total area was more than 10−10, then the polygon Ei-F1-F2-,…,-FK-Ei + 1 was taken as a background cell to increase the accuracy of integration. The nodes {F1, F2,…,FK} were called the background integral points. As shown in Figure 6, the triangles △A1-D1-C1, △C1-D2-C2, △C2-D3-C3, △C3-D4-C4, △C4-D5-C5, △C5-D6-A2, △A2-D7-C6, △C6-D8-C7, △C7-D9-C8, △C8-D10-C9, △C9-D11-C10, and △C10-D12-A1 were background cells. The nodes D1,D2,…,D12 were background integral points. The folded lines A1-D1-C1-D2-C2-D3-C3-D4-C4-D5-C5-D6-A2-D7-C6-D8-C7-D9-C8-D10-C9-D11-C10-D12-A1 were the revised approaching boundary. Through the above steps, the initial approaching boundary and integration domains were first revised.After that, the nodesA1,A2,…,An, C1,C2,…,Cl, and D1,…,Dk in each irregular element were figured out, and they were sorted from small to large in the view of the x coordinates, recorded as {G1, G2,…,Gll}. Every two adjacent nodes, Gj (xj, yj), Gj + 1 (xj + 1, yj + 1), were connected into a straight line Gj-Gj + 1, and the midpoint of straight line Gj-Gj+1 was recorded as Hj. A node Pj (xj′, yj′) that was on the real boundary and close to Hj needed to be figured out, in which yi′ and xi′ satisfied the following equations:(37)yi′=yi+yi+12,fxi′,yi′=0,where f is the function of the real boundary. The nodes Pi (xi′, yi′), Gi (xi, yi), and Gi + 1 (xi + 1, yi + 1) were connected to construct a triangle, △Gi-Gi + 1-Pi. If the ratio of the area of triangle △Gi-Gi + 1-Pi to the total area was more than 10−10, then the triangle △Gi-Gi + 1-Pi was taken as a background cell to increase the accuracy of integration. The node Pi (xi′, yi′) was called the background integral point, which was not used for interpolation and would not increase the degrees of freedom. In Figure 7, the nodes D9, D10, P1, and P2 are background integral points, while the triangles △B12-D9-C8, △C8-D10-B13, △D9-P1-C8, and △C8-P2-D10 are background cells.Figure 7
Background integral points and background cells.The above process was repeated to seek out all background integral points and background cells. In this way, a good approximation for the physical boundaries and accurate integration domains can be regained. Finally, the integration of all background cells was integrated into the global stiffness matrixK and vector f. Therefore, good mesh integration accuracy was achieved without increasing the degrees of freedom as the background integral points were not used to interpolate.


## 6. Numerical Examples

### 6.1. Clamped Circular Plate under Distributed Load

A clamped circular plate under a distributed load of intensityp is shown in Figure 8. The analytical solutions of the deflection and moment for this numerical example [39] can be written as follows:(38)w=px2+y2−r264D,where D = 1, p = 1, μ = 0.3, and r = 1.Figure 8
A clamped circular plate.The deflection error is defined as follows:(39)werr=w¯−ww×100%,where w and w¯ are the analytical and numerical solutions of the deflection, respectively.The meshes used for the proposed NMM and the traditional NMM are shown in Figure9. The meshes used for various FEMs are shown in Figure 10.Figure 9
Meshes with different element sizes of a clamped circular plate divided for the proposed NMM and the traditional NMM; from left to right: 0.5 × 0.5, 0.3 × 0.3, 0.25 × 0.25, 0.2 × 0.2, and 0.15 × 0.15.Figure 10
Meshes with different elements of a clamped circular plate divided for various FEMs; from left to right: 48, 108, 192, 432, and 768.Many commercial software programs developed based on the FEM have been successfully applied to various fields, and the theory of simulation for plate and shell problems is also well established. To illustrate the efficiency of the proposed NMM developed in this study, the element of type Shell63 in ANSYS 10.0 developed by ANSYS Inc. (Pittsburgh, PA, USA) was selected. The element of type Shell63 is a plate element with six parameters: the displacements in thex, y, and z axes and the rotations about these three axes. Because the distributed load is perpendicular to the neutral surface, the displacements in the x and y axes and the rotation about the z axis can be considered to be equal to 0. Thus, the element type Shell63 can simulate thin plate bending problems.We used the element of type Shell63 in ANSYS 10.0, the proposed NMM, the traditional NMM, and the analytical method to analyze this circular plate. And the results are reported in Table1. In the first and fourth columns of Table 1, the length represents the largest element diagonal. Figure 11 shows the deflections achieved with the proposed NMM and the analytical method. The results achieved with the proposed NMM compared well with the analytical solutions. The convergence processes of the central deflection of the clamped circular plate, under distributed load, calculated with the proposed NMM, were oscillatory. Figure 12 shows the central deflections (left) and error comparisons (right) calculated with different methods. It can be found that the deflections achieved with the proposed NMM were much better than those calculated with the traditional NMM, the FEM, RPAQ [40, 41], and TACQ [42].Table 1
Central deflections of a clamped circular plate under distributed load calculated with different methods.
Length (10−1m)/(element number)Central deflections (10−2 m)Length (10−1m)/(element number)Central deflections (10−2 m)The traditional NMMThe proposed NMMFEMRPAQ [40, 41]TACQ [42]7.0711/(16)1.6530 (5.79%)1.5997 (2.38%)7.368/(48)1.5819 (1.24%)1.264 (19.1%)—4.2426/(52)1.5717 (0.59%)1.5541 (0.53%)4.3051/(108)1.5715 (0.58%)1.484 (5.02%)1.4790 (5.34%)3.5355/(60)1.5629 (0.03%)1.5614 (0.07%)2.3301/(192)1.5672 (0.3%)1.543 (1.25%)1.5422 (1.30%)2.8284/(88)1.5553 (0.46%)1.5627 (0.01%)1.5696/(432)1.5644 (0.12%)——2.1213/(164)1.5628 (0.02%)1.5625 (0%)1.1814/(768)1.5635 (0.06%)——Analytical solution1.5625Analytical solution1.5625Note. “—” indicates that no specific results are given.Figure 11
Deflection of a clamped circular plate under a distributed load calculated with the proposed NMM with 0.15 × 0.15 mesh and the analytical method.Figure 12
Central deflections (left) and error comparisons (right) of a clamped circular plate under distributed load calculated with different methods.


### 6.2. Simply Supported Annular Plate under Distributed Load

An annular plate depicted in Figure13(a) under the distributed load of intensity p was studied here, with its outer edge being simply supported and its inner edge being free. The analytical solution of deflection and moment for this numerical example [39] can be expressed as follows:(40)w=pa464D−1−x2+y2a22+2a11+u1−x2+y2a2−4a2β21−uInx2+y2a2,where p=0.001, μ = 0.3, and a1=3+μ3−β241+μβ2κ,a2=3+μ+41+μκ,κ=β2/1−β2Inβ, β = b/a, with a = 5 and b = 1 as the outer and inner radii of the annular plate, respectively. The meshes used for the traditional NMM and the proposed NMM are shown in Figure 13(b), and the meshes used for the FEM are shown in Figure 13(c).Figure 13
An annular plate with its outer edge being simply supported and its inner edge being free (a); 0.2 × 0.2 mesh of the plate divided for the proposed NMM and the traditional NMM (b); the mesh with 3200 elements of the plate divided for the FEM (c).
(a)(b)(c)We used the element of type Shell63 in ANSYS 10.0, the proposed NMM, the traditional NMM, and the analytical method to analyze the simply supported annular plate. Figure14 displays the deflections (left) and deflection errors (right) achieved with the proposed NMM with 0.2 × 0.2 mesh. Figure 15 presents the deflections on the line y = 0 achieved with the proposed NMM, the traditional NMM, the element of type Shell63 in ANSYS 10.0, and the analytical method. It can be found that the deflections calculated with the proposed NMM were compared well with the analytical solutions. Additionally, the deflections achieved with the proposed NMM were slightly better than those calculated with the traditional NMM and the FEM.Figure 14
Deflections (a) and deflection errors (b) of a simply supported annular plate under a distributed load calculated with the proposed NMM with 0.2 × 0.2 mesh.
(a)(b)Figure 15
Deflection comparisons of a simply supported annular plate under a distributed load calculated with different methods.Compared with the results of the above examples, the validity of the proposed NMM used to solve the convergence and numerical accuracy of ACM element with irregular physical boundaries was verified. Additionally, the results calculated with the proposed NMM introduced in this study were more accurate than those achieved with the traditional NMM and the FEM discussed in the examples, since the proposed NMM incorporated the regular meshes, the background cells, and the mixed-order local displacement function in the bending analysis of nonconforming thin plates. The regular meshes were selected to improve convergence performance. The background cells were employed to obtain good integration accuracy. The mixed-order local displacement function was taken to improve the interpolation accuracy. In this way, the best numerical accuracy was achieved, and the convergence was correspondingly reached.


## 6.1. Clamped Circular Plate under Distributed Load

A clamped circular plate under a distributed load of intensityp is shown in Figure 8. The analytical solutions of the deflection and moment for this numerical example [39] can be written as follows:(38)w=px2+y2−r264D,where D = 1, p = 1, μ = 0.3, and r = 1.Figure 8
A clamped circular plate.The deflection error is defined as follows:(39)werr=w¯−ww×100%,where w and w¯ are the analytical and numerical solutions of the deflection, respectively.The meshes used for the proposed NMM and the traditional NMM are shown in Figure9. The meshes used for various FEMs are shown in Figure 10.Figure 9
Meshes with different element sizes of a clamped circular plate divided for the proposed NMM and the traditional NMM; from left to right: 0.5 × 0.5, 0.3 × 0.3, 0.25 × 0.25, 0.2 × 0.2, and 0.15 × 0.15.Figure 10
Meshes with different elements of a clamped circular plate divided for various FEMs; from left to right: 48, 108, 192, 432, and 768.Many commercial software programs developed based on the FEM have been successfully applied to various fields, and the theory of simulation for plate and shell problems is also well established. To illustrate the efficiency of the proposed NMM developed in this study, the element of type Shell63 in ANSYS 10.0 developed by ANSYS Inc. (Pittsburgh, PA, USA) was selected. The element of type Shell63 is a plate element with six parameters: the displacements in thex, y, and z axes and the rotations about these three axes. Because the distributed load is perpendicular to the neutral surface, the displacements in the x and y axes and the rotation about the z axis can be considered to be equal to 0. Thus, the element type Shell63 can simulate thin plate bending problems.We used the element of type Shell63 in ANSYS 10.0, the proposed NMM, the traditional NMM, and the analytical method to analyze this circular plate. And the results are reported in Table1. In the first and fourth columns of Table 1, the length represents the largest element diagonal. Figure 11 shows the deflections achieved with the proposed NMM and the analytical method. The results achieved with the proposed NMM compared well with the analytical solutions. The convergence processes of the central deflection of the clamped circular plate, under distributed load, calculated with the proposed NMM, were oscillatory. Figure 12 shows the central deflections (left) and error comparisons (right) calculated with different methods. It can be found that the deflections achieved with the proposed NMM were much better than those calculated with the traditional NMM, the FEM, RPAQ [40, 41], and TACQ [42].Table 1
Central deflections of a clamped circular plate under distributed load calculated with different methods.
Length (10−1m)/(element number)Central deflections (10−2 m)Length (10−1m)/(element number)Central deflections (10−2 m)The traditional NMMThe proposed NMMFEMRPAQ [40, 41]TACQ [42]7.0711/(16)1.6530 (5.79%)1.5997 (2.38%)7.368/(48)1.5819 (1.24%)1.264 (19.1%)—4.2426/(52)1.5717 (0.59%)1.5541 (0.53%)4.3051/(108)1.5715 (0.58%)1.484 (5.02%)1.4790 (5.34%)3.5355/(60)1.5629 (0.03%)1.5614 (0.07%)2.3301/(192)1.5672 (0.3%)1.543 (1.25%)1.5422 (1.30%)2.8284/(88)1.5553 (0.46%)1.5627 (0.01%)1.5696/(432)1.5644 (0.12%)——2.1213/(164)1.5628 (0.02%)1.5625 (0%)1.1814/(768)1.5635 (0.06%)——Analytical solution1.5625Analytical solution1.5625Note. “—” indicates that no specific results are given.Figure 11
Deflection of a clamped circular plate under a distributed load calculated with the proposed NMM with 0.15 × 0.15 mesh and the analytical method.Figure 12
Central deflections (left) and error comparisons (right) of a clamped circular plate under distributed load calculated with different methods.


## 6.2. Simply Supported Annular Plate under Distributed Load

An annular plate depicted in Figure13(a) under the distributed load of intensity p was studied here, with its outer edge being simply supported and its inner edge being free. The analytical solution of deflection and moment for this numerical example [39] can be expressed as follows:(40)w=pa464D−1−x2+y2a22+2a11+u1−x2+y2a2−4a2β21−uInx2+y2a2,where p=0.001, μ = 0.3, and a1=3+μ3−β241+μβ2κ,a2=3+μ+41+μκ,κ=β2/1−β2Inβ, β = b/a, with a = 5 and b = 1 as the outer and inner radii of the annular plate, respectively. The meshes used for the traditional NMM and the proposed NMM are shown in Figure 13(b), and the meshes used for the FEM are shown in Figure 13(c).Figure 13
An annular plate with its outer edge being simply supported and its inner edge being free (a); 0.2 × 0.2 mesh of the plate divided for the proposed NMM and the traditional NMM (b); the mesh with 3200 elements of the plate divided for the FEM (c).
(a)(b)(c)We used the element of type Shell63 in ANSYS 10.0, the proposed NMM, the traditional NMM, and the analytical method to analyze the simply supported annular plate. Figure14 displays the deflections (left) and deflection errors (right) achieved with the proposed NMM with 0.2 × 0.2 mesh. Figure 15 presents the deflections on the line y = 0 achieved with the proposed NMM, the traditional NMM, the element of type Shell63 in ANSYS 10.0, and the analytical method. It can be found that the deflections calculated with the proposed NMM were compared well with the analytical solutions. Additionally, the deflections achieved with the proposed NMM were slightly better than those calculated with the traditional NMM and the FEM.Figure 14
Deflections (a) and deflection errors (b) of a simply supported annular plate under a distributed load calculated with the proposed NMM with 0.2 × 0.2 mesh.
(a)(b)Figure 15
Deflection comparisons of a simply supported annular plate under a distributed load calculated with different methods.Compared with the results of the above examples, the validity of the proposed NMM used to solve the convergence and numerical accuracy of ACM element with irregular physical boundaries was verified. Additionally, the results calculated with the proposed NMM introduced in this study were more accurate than those achieved with the traditional NMM and the FEM discussed in the examples, since the proposed NMM incorporated the regular meshes, the background cells, and the mixed-order local displacement function in the bending analysis of nonconforming thin plates. The regular meshes were selected to improve convergence performance. The background cells were employed to obtain good integration accuracy. The mixed-order local displacement function was taken to improve the interpolation accuracy. In this way, the best numerical accuracy was achieved, and the convergence was correspondingly reached.


## 7. Conclusions

After the comparisons between the characteristics of several conforming and nonconforming elements were conducted in detail, the ACM element was adopted in this paper. The biggest issue in current analytical methods is that the convergence, good numerical accuracy, and high computing efficiency of ACM element with irregular physical boundaries cannot be achieved simultaneously. In order to solve the problem, the mixed-order NMM with background cells for integration was developed and introduced in this paper. The discussed examples showed that the numerical solutions achieved with the proposed NMM rapidly converged to the analytical solutions and their accuracy was vastly superior to that achieved with the FEM and the traditional NMM, which verified the proposed NMM.



---

*Source: 6681214-2020-12-12.xml*