A method for treating nonconforming thin plate bending problems

Finite Elements in Analysis and Design 47 (2011) 200–207

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Review

A method for treating nonconforming thin plate bending problems Zhongxu Tian a,n, Delin Wang b a b

Key Laboratory for Power Machinery and Engineering of Ministry of Education, Shanghai Jiaotong University, Shanghai 200240, PR China Department of Aerospace Engineering, Texas A&M University, College Station, TX 77843, USA

a r t i c l e in f o

a b s t r a c t

Article history: Received 12 June 2009 Received in revised form 3 September 2010 Accepted 6 September 2010

Weak form for the thin plate bending problems is derived and the boundary parameters included. The weak form is useful in dealing with and understanding about nonconforming problems in finite element methods and other numerical methods. It can be observed that finite element method (FEM) can also be formulated from the same weak form. In addition, a new discretization scheme is introduced, which is proven by numerical results to be highly accurate and has no nonconforming problems. Crown Copyright & 2010 Published by Elsevier B.V. All rights reserved.

Keywords: Weak form Thin plate bending Boundary conditions Nonconforming Discretization scheme Finite element method

Contents 1. 2. 3. 4. 5. 6.

7.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 Weak forms for thin plate bending problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 Treatment of boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Finite element method and nonconforming problems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Another discrete schemes to avoid nonconforming problems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 6.1. Patch test for Triangle Element 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 6.2. Patch test for triangle element 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 6.3. Patch test for quadrilateral element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 6.4. Evaluation of displacement for rectangular plates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207

1. Introduction A plate is a flat structural element with planform dimensions much larger than its thickness and is subjected to loads that cause bending deformation in addition to stretching. Up to now, the plate elements can be categorized into three types: (1) the C1 continuity elements based on Poisson–Kirchhoff plate theory, (2) the C0 continuity elements based on Reissner–Mindlin plate theory and (3) the hybrid elements [1]. The Poisson–Kirchhoff theory of plates was the first displacement-based finite element approach to be used in the bending analysis [2–5]. The theory

n

Corresponding author. E-mail address: [email protected] (Z.X. Tian).

does not include the shear deformation and is valid only for thin plates. Moreover, the theory requires C1 continuity to secure the interelement compatibility which would include a difficult task of maintaining the continuity of the slope of displacement across the element boundaries [6]. The difficulty of imposing C1 continuity has resulted in many alternative approaches to the problem. Two basic alternatives are as follows. First, the Reissner–Mindlin plate theory which requires only C0 continuity and accommodates transverse shear strains was introduced in the formulation of plate bending. The theory works well for thick plates. However, as the plate thickness is decreased, the shear terms dominate in the stiffness matrix and the so-called ‘‘locking’’ phenomenon occurs. In the past decades much attention has been forcused on the investigation of shear locking and relevant numerical difficulties, which can be found in

0168-874X/$ - see front matter Crown Copyright & 2010 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2010.09.002

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

Refs. [7–12]. A second alternative to overcome the C1 continuity problem is to modify the potential energy by adding Lagrangian multiplier terms and thereby introducing independent field parameters subject to variation in the variational treatment. The basic idea behind this concept is to include the strain energy produced on the interelement boundaries from normal slope discontinuities in the formulation of the total energy balance of the structure. Elements of this kind which is called hybrid elements have been suggested in [13–18] among others. In recent decades, more and more people studied FEM using weak form equations, and FEM can be seen as an weighted residual-Galerkin method. For weak form equations, the shape functions and test functions must be carefully selected so that the integrations must be valued [19]. Furthermore, finite volume method (FVM) also using weak forms, and seldom using interpolation functions, which is commonly used in the computational fluid mechanics [20–23]. However, in computational solid mechanics FVM is not widely used and the relationship between FVM and FEM remains unclear. Onate et al. [24] claimed that FVM and FEM share in common the concepts such as mesh discretization and interpolation across elements. As an extension of traditional FVM, cell centered finite volume method (CFVM) needs no special treatment of the interelement C1 continuity conditions [25]. Fallah [26] developed a new form of finite volume method based on Mindlin–Reissner plate theory for the analysis of transversely loaded plates which appears to be locking free in the deformation prediction. A novel finite volume based discretization method was presented in [27] for loaded two dimensional structures, which incorporates rotation variables in addition to the displacement variables in the finite volume based structural analysis procedures. While giving FVM schemes of thin plate bending problems, there are still some problems kept unclear. Firstly, in order to struggle with four-order difference operators, we need to set up interpolation functions to express higher order derivatives of deflection, whether the interpolation should still be C1. Furthermore, how to express the force boundary conditions which are given with determined bending moments and transverse shearing forces. In this paper, a new method like FVM is introduced to solve thin plate bending problems based on weak forms of difference equations. The relationship between FEM and the new method, and nonconforming problems are investigated through some instances. Furthermore, boundary conditions are treated in a way which is easily applicable in engineering. The whole paper is structured as follows: in the first section we outline the weak form for thin plate bending problem. In the next three sections we address the treatment of boundary conditions and the singularity problems. Section 6 provides the illustration of the discrete scheme and relevant numerical results. The main observations are summarized in the concluding remarks. 2. Weak forms for thin plate bending problem The governing equation of plate bending problem is ! @2 Mxy @2 M @2 Mx D þ þ 2 q ¼ 0 8ðx, yÞ A O @x@y @y2 @x2

ð1Þ

where D is constant relating to rigidity coefficient E, Poisson’s ratio v and the thickness t of the plates. Mx, Mxy and My are the moments per unit length which can be put in terms of transverse deflection w as follows: Mx ¼ D

! @2 w @2 w þv 2 , 2 @x @y

My ¼ D

! @2 w @2 w þv 2 , 2 @y @x

Mxy ¼ Dð1vÞ

@2 w @x@y

ð2Þ

201

Assuming that a test function P is applied on an arbitrary domain O0 C O, and zero outside of domain O0 . We can obtain the weak form to Eq. (1): ! ZZ ZZ @2 Mx @2 Mxy @2 My P dxdy þ  D þ2 qP dx dy ¼ 0 ð3Þ 2 2 @x@y @x @y O0 O0 where w is usually called the sample function in the moments as shown in (2) while P the test function. By applying Green’s formula twice on Eq. (3), we obtain  Z  Z @P @P þMns ds Qn V ds þ Mn  @n @s @O0 @O0 ! ZZ ZZ @2 P @2 P @2 P þMy 2 dx dy  Mx 2 þ2Mxy qP dx dy ¼ 0 ð4Þ @x@y @x @y O0 @O0 where qO0 denotes the boundary of O0 , and Mn, Mns and Qn are bending moment, twisting moment and transverse shearing force, respectively, which can be expressed in the following forms: Mn ¼ n2x Mx þn2y My þ 2nx ny Mxy Mns ¼ nx ny Mx þ nx ny My þ ðn2x n2y ÞMxy Qn ¼ Qx nx þ Qy ny

ð5Þ

where nx, ny are cosine directions of outward normal of qO0 . Transverse shearing forces Qx and Qy are related to the bending and twisting moments Mx, Mxy and My in the following way: Qx ¼

@Mx @Mxy þ , @x @y

Qy ¼

@Mxy @My þ @x @y

In the end, Eq. (3) can be written as Z Z @P @P A2 dx dy A3 dx dy A1 P ds @O0 @Ou @x @Ou @y ZZ ZZ 1 þ A dx dy qP dx dy ¼ 0 D Ou Ou

ð6Þ

Z

ð7Þ

where A1 ¼

A2 ¼

A3 ¼

A¼

@3 w @3 w @3 w @3 w nx þ nx þ 3 ny þ 2 ny @x3 @x@y2 @y @x @y ! @2 w @2 w @2 w ny þ v 2 nx þ ð1vÞ @x@y @x2 @y ! @2 w @2 w @2 w nx þ v ny þ ð1vÞ 2 2 @x@y @y @x

ð8Þ

ð9Þ

ð10Þ

@2 P @2 w @2 P @2 w @2 P @2 w @2 P @2 w @2 P @2 w þ 2 2 þv 2 2 þ v 2 2 þ 2ð1vÞ @x@y @x@y @x2 @x2 @y @y @x @y @y @x ð11Þ

Eq. (7) is the one employed in discretization scheme. Patch functions are usually used in order to fit irregular domains. In Eq. (7), expression A includes second order derivatives to coordinates, which goes to infinity when the shape functions and test functions are not C1 in domain O0 , and its integration are incapable of being evaluated [19]. Two methods can be the alternatives to overcome the problem: (1) Appropriate C1 type sample and test function are chosen to make sure that A is finite and integrable in Eq. (7). Conforming element in FEM is an example of such method. (2) Test function P is chosen to be smooth and linear over domain O0 so that A is zero on O0 . This is the method used in the current work as shown in later section.

202

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

3. Treatment of boundary conditions ZZ In further proceeding with the above weak form (4), it will be convenient to use Mx, Mxy and My as internal variables while dealing with the boundary conditions. If O0 touches qO (Fig. 1), qO0 should be divided into two parts qO0 \qO and qO0  qO0 \qO, then Eq. (4) becomes   Z Z @P @P þ Mns ds Qn V ds þ Mn  @n @s @O0 \@O @O0 \@O ! ZZ ZZ @2 P @2 P @2 P  þ My 2 dx dy Mx 2 þ 2Mxy qP dx dy 0 @x@y @x @y O @O0   Z Z @P @P þ Mns ds ¼ 0 ð12Þ  Q n V ds þ Mn 0 0 @n @s @O@O \@O @O@O \@O

where Q n , M ns , M nn are forces acting on qO0 \qO. It can be noted that boundary conditions can be taken care of easily according to (12), and the disposal method is similar to that of FEM.

Replacing P in Eq. (7) with N, we have by Eq. (13): ZZ 1 A dx dy ¼ qN dxdy D Ou Ou

ð14Þ

where A¼

@2 N @2 w @2 N @2 w @2 N @2 w @2 N @2 w @2 N @2 w þ 2 þv 2 þv 2 þ 2ð1vÞ 2 2 2 2 2 @x@y @x@y @x @x @y @y @x @y @y @x ð15Þ

Since slopes of N and w are C1, A is finite in O0 and Eq. (14) can be written as 5 ZZ 5 ZZ X 1X A dx dy ¼ qN dxdy ð16Þ D i ¼ 1 Oi Oi i¼1 It can be noted that the implementation of FEM is based on Eq. (16). The summation in the left hand side of the equation reflects in FEM the assembling of element stiffness matrix into general stiffness matrix. Eq. (16) is the one employed in the formulation of conforming elements in FEM. However, for nonconforming elements, the use of Eq. (16) might result in errors due to the following two reasons:

4. Finite element method and nonconforming problems Assume that the deflection w is interpolated by patch functions which related to some finite element mesh and nodal parameters, nodal deflections and slopes. For a typical patch as shown in Fig. 2, an arbitrary node O is the vertex shared by elements Oi (i¼1,y,5), which consist of domain O0 . Let N denote the base function corresponding to a nodal parameter of O, named shape function in FEM. Suppose deflection function is conforming, and deflection w should be C1 across element boundaries for thin bending problems. Being one base function of w, N must be a C1 function. As w is C1 at the interelement boundaries, the deflections and slopes on the qO0 should not be affected by nodal parameters of node O, so N should satisfies N¼

@N @N ¼ ¼0 @x @y

ðon @OuÞ

ð13Þ

∂Ω Ω

Ω′ ∂Ω′

(a) The value of A goes to infinity [19] at the interelement boundary and Eq. (14) cannot be substituted with Eq. (16). (b) Without condition (13), integrates along boundary qO0 in Eq. (7) are not zero and Eq. (14) do not exist.

5. Another discrete schemes to avoid nonconforming problems To circumvent the problem, a new discrete scheme is in introduced as follows. The test function P can be chosen to be smooth and only linear in O0 so that the integral of A in O0 vanishes automatically, and nonconforming problems no longer exist. As it is illustrated in Fig. 3 for an arbitrary node O, the test functions are chosen to be 1, x xi and y yi, respectively, and (xi, yi) are coordinates of node O. O0 is chosen in a way that qO0 pass through the centroids of all triangle or quadrilateral elements sharing node O as vertex and the midpoints of the element edges. Thus by substituting P¼1, x xi and y yi into (7), respectively, we have the following discrete form: RR 8 R A dxdy ¼ Ou qdx dy D > < R@Ou 1 R RR D @Ou ðxxi ÞA1 dxdyD @Ou A2 dxdy ¼ Ou qðxxi Þ dx dy ð17Þ R RR R > :D ðyy ÞA dxdyD A dxdy ¼ qðyy Þ dx dy 1 3 i i @Ou @Ou Ou The numerical solution can be obtained after all the discrete equations on every node in domain O are set up. The same process

Fig. 1. Domain O and domain O0 .

4

Ω4

5

3

Ω3 Ω2

O Ω1

Ω5 1

Fig. 2. Domain O0 .

O P2

2 P1

Fig. 3. Domain O0 for the new discrete scheme.

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

as FEM can introduced to implement the new scheme, except the stiffness matrix are different. The new stiffness matrix can be computed according to left side of Eq. (17) and by borrowing Gauss integration method. In Eq. (17), only boundary integrations exist, and qO0 pass through the interior of the concerned elements, so no infinite items needed to be evalued and nonconforming interpolation functions might be workable in Eq. (17). Discrete form (17) has another type of expression as follows if P¼1, x xi and y yi are substituted into Eq. (4), respectively 8 R RR  @O0 Qn ds @O0 q dx dy ¼ 0 > > >   > > R > <  0 Q ðxx Þ dsþ R 0 M @ðxxi Þ þ M @ðxxi Þ dsRR 0 qðxx Þ dx dy ¼ 0 n ns i i @O n @O @O @n @s   > > R R R R > @ðyy Þ @ðyy Þ > i i > ds @O0 qðyyi Þ dx dy ¼ 0 > :  @O0 Qn ðyyi Þ dsþ @O0 Mn @n þ Mns @s

203

It can be noted that the three equations of (18) stands for the force equilibrium conditions in z direction and the moment equilibrium around the axes in y and x directions passing through node i. Being another form of Eq. (18), Eq. (17) are the equilibrium conditions for domain O0 , so the physical meaning of the discrete scheme given in Eq. (17) has more clarity than that of FEM.

6. Numerical examples In this section, several test problems are analyzed to ascertain the accuracy and convergence of the new discrete scheme as reflected in Eq. (17) or (18). 6.1. Patch test for Triangle Element 1

ð18Þ Here, Triangle Element 1 is given by the new discrete scheme Eq. (18), and the element interpolation function is chosen to be

and, we have @ðxxi Þ @ðxxi Þ þ Mns ¼ Mn nx Mns ny @n @s @ðyyi Þ @ðyyi Þ þ Mns ¼ Mn ny þ Mns ny Mn @n @s

w ¼ a1 þ a2 x þ a3 y þ a4 x2 þ a5 xy þ a6 y2 þ a7 x3 þ a8 y3 þ a9 ðx2 y þxy2 Þ

Mn

ð20Þ ð19Þ In this example triangle element 1 is considered for the patch test in comparison with finite element sharing the same element interpolation function. The constants in Eq. (20) can be determined by substituting the nodal parameters and coordinates into the equation and thus the interpolation function of displacement w can be defined. It needs to be pointed out that the function w is only C0 type which can be verified by investigating the deflections and slopes on common boundaries between any two neighbored elements. An arbitrary mesh consisting of triangle elements is shown in Fig. 4 with nodal coordinates tabulated in Table 1. The boundary deflections and slopes are given to be an arbitrary function

8

7

6

5

4

3

w ¼ 0:17 þ0:02x þ 0:04y þ0:05x2 þ0:08xy þ 0:1y2 þ0:2x3 þ 0:22y3 þ 0:18ðx2 y þxy2 Þ

ð21Þ

2

1 Fig. 4. Triangular-element mesh for patch test.

Table 1 Nodal coordinates of triangle elements. Node

1

2

3

4

5

6

7

8

X Y

0.00 0.00

0.50 0.00

0.20 0.18

0.38 0.20

0.35 0.35

0.23 0.37

0.00 0.50

0.50 0.50

Function (21) follows the form of element interpolation function (20) and the coefficients are given arbitrarily. The numerical results of deflections and slopes of internal nodes are listed in Table 2 in comparison with the analytical solutions and FEM results using the same interpolation function expressed by (20). It can be seen from Table 2 that even for nonconforming element interpolation function, the triangle element based on the present discrete scheme can pass the patch test which is much stricter than that usually used in FEM. It

Table 2 Numerical results of patch test for triangle element 1. Node

Node 3

Node 4

Deflections and slopes

w

@w @x

@w @y

w

@w @x

@w @y

Triangle element 1 FEM Exact

 0.1465336  0.1443589  0.1465336

0.0971920 0.0904111 0.097192

0.1335440 0.1267450 0.133544

 0.1164312  0.1135488  0.1164312

0.1952000 0.2059242 0.1952

0.1901520 0.2074295 0.190152

Node

Node 5

Deflections and slopes

w

@w @x

@w @y

w

@w @x

@w @y

Triangle element 1 FEM Exact

 0.0873825  0.0841229  0.0873825

0.2226500 0.2332224 0.22265

0.2850000 0.2756933 0.285

 0.1046891  0.1041765  0.1046891

0.1596180 0.1737998 0.159618

0.2629120 0.2615619 0.262912

Node 6

204

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

should be noted that the constants in Eq. (20) cannot be determined for some right triangle, so element 1 is not fit for practical application.

Following the form of Eq. (23), we specify the deflections and slopes for nodes on the plate boundary according to the following function:

6.2. Patch test for triangle element 2

w ¼ 0:17 þ 0:02x þ 0:04yþ 0:05x2 þ0:08xy þ 0:1y2 þ 0:2x3 þ 0:18x2 y þ 0:14xy2 þ 0:22y3

The triangle element 2 is also pushed forward using the new discrete scheme Eq. (18) and the interpolation function is selected the same way as that of BCIZ element well known in FEM field [5]. It is very easy to verify as shown in the following numerical tests that triangular element 2 can pass the patch test which BCIZ element cannot. Numerical results are listed in Table 3 of the same patch test as the boundary deflections and slopes are given to be an arbitrary function w ¼ 0:17 þ0:02x þ 0:04y þ0:05x2 þ 0:08xy þ0:1y2

ð26Þ

Function (26) can be expressed by (23) and the coefficients are selected arbitrarily. The computed deflections and slopes for the internal nodes in comparison with the finite element sharing the same element interpolation function are listed in Table 4. It is obvious that the quadrilateral element based on the new discrete scheme can pass the patch test. Moreover, the new discrete scheme causes no loss in accuracy even if there is a nonconforming problem.

ð22Þ

8

7 6.3. Patch test for quadrilateral element

6

The test which was introduced by Robinson [28] has an assembly of five quadrilateral elements as shown in Fig. 5. The element interpolation function of quadrangle element is selected in the following form for the new discrete scheme given in Eq. (17)

5

4

3

w ¼ a1 þ a2 x þ a3 y þ a4 x2 þ a5 xy þ a6 y2 þ a7 x3 þ a8 x2 y þ a9 xy2 þ a10 y3 þ a11 xy3 þ a12 x3 y

ð23Þ

1

2

Fig. 5. Quadrilateral element mesh for patch test.

Fig. 6 shows two neighbored elements, a and b, shared boundary AB which can be expressed with y ¼ k1 x þ k2

x

ð24Þ

According to element function (23), the deflections at boundary AB can be given in the following forms: waAB ¼ b1 þb2 x þ b3 x2 þ a4 x3 þ a5 x4

B

wbAB ¼ c1 þ c2 x þc3 x2 þc4 x3 þc5 x4

ð25Þ

a

where waAB and wbAB are the deflections on AB determined with deflection functions of elements a and b, respectively. Along direction s, there are only four common parameters for the two elements, wA, wB, (qw/qt)A and (qw/qt)B, to constrain the five parameters of waAB and wbAB , so waAB and wbAB cannot be kept consistent, and the quadrangle element is neither C0 nor C1.

b

τ

A y Fig. 6. Two neighbored quadrilateral elements.

Table 3 Numerical results of patch test for triangle element 2. Node

Node 3

Node 4

Deflection and slopes

w

@w @x

@w @y

w

@w @x

@w @y

Triangle element 2 BCIZ element Exact

 0.1506800  0.1516317  0.15068

0.0544000 0.0514140 0.0544

0.0920000 0.0913744 0.092

 0.1371000  0.1377298  0.1371

0.0740000 0.0760691 0.074

0.1104000 0.1113600 0.1104

Node

Node 5

Deflection and slopes

w

@w @x

@w @y

w

@w @x

@w @y

Triangle element 2 BCIZ element Exact

 0.1208250  0.1214446  0.120825

0.0830000 0.0859902 0.083

0.1380000 0.1400165 0.138

 0.1274570  0.1281440  0.127457

0.0726000 0.0724478 0.0726

0.1324000 0.1343920 0.1324

Node 6

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

205

Table 4 Numerical results of patch test for quadrilateral element. Node

Node 3

Node 4

Deflections and slopes

w

@w @x

@w @y

w

@w @x

@w @y

Present element FEM Exact

 0.1455938  0.14623  0.1455938

0.0958960 0.096529 0.095896

0.1306640 0.13303 0.130664

 0.1170392  0.11677  0.1170392

0.193600 0.20195 0.1936

0.1840720 0.18975 0.184072

Node

Node 5

Deflection and slopes

w

@w @x

@w @y

w

@w @x

@w @y

Present element FEM Exact

 0.0890975  0.088561  0.0890975

0.2177500 0.21826 0.2177500

0.2752000 0.27630 0.2752000

 0.1059486  0.10555  0.1059486

0.1541420 0.15879 0.1541420

0.2561040 0.26441 0.2561040

Node 6

C

C

N=2

N=2

N=2

N=2

Mesh A

Mesh B

Fig. 7. Two types of triangular element meshes defined on the quarter-plate.

Table 5 The numerical solution wC D=FL2 from mesh type A versus the exact one where L denotes the length of the plate and F the concentrated load. Boundary

Simply supported

Clamped

Mesh

22

44

88

16  16

22

44

88

16  16

Element 2 BCIZ element Veubeke element exact

0.01153 0.01231 0.01062

0.01155 0.01205 0.01127

0.01158 0.01199 0.01150

0.01159 0.01198 0.01157

0.005248 0.005837 0.004302

0.005522 0.005825 0.005249

0.005579 0.005799 0.005499

0.005602 0.005792 0.005579

0.01160

0.005612

Table 6 The numerical solution wC D=FL2 from mesh type B versus the exact one where L denotes the length of the plate and F the concentrated load. Boundary

Simply supported

Clamped

Mesh

22

44

88

16  16

22

44

88

16  16

Element 2 BCIZ element Veubeke element Exact

0.01225 0.01375 0.01048

0.01186 0.01244 0.01123

0.01170 0.01194 0.01149

0.01164 0.01174 0.01157

0.005509 0.006564 0.004215

0.005681 0.006136 0.005159

0.005658 0.005837 0.005476

0.005633 0.005706 0.005573

0.01160

6.4. Evaluation of displacement for rectangular plates Square clamped and simply supported plates subjected to a concentrated central load are analyzed using triangular element 2 in comparison with the well known BCIZ element [5] and Veubeke [3]. Two types of triangle element meshes, mesh A and mesh B as

0.005612

illustrated in Fig. 7, are considered on a quarter-plate. The mesh size is defined by the number of divisions per mid-sides on the quarter-plate which for the current analysis are 2  2, 4  4, 8  8, and 16  16, respectively. The transverse displacement w at point C is reported in Tables 5 and 6. For both the simply supported and clamped square plates under a central point load, the relative

206

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

9

element 2 BCIZ Element Veubeke element

8

6

Relative erroe %

Relative error %

7

5 4 3 2 1 0 0

2

4

6

8

10

12

14

16

18

26 24 22 20 18 16 14 12 10 8 6 4 2 0 -2

Element 2 BCIZ Element Veubeke Element

0

2

25

Relative error %

8

10

12

14

16

18

Fig. 11. Relative error Z in transverse displacement at point C versus mesh size for clamped square plate and mesh type B.

error in central displacements defined as  C  w wC  EX   Z ¼  C   100% wEX

Element 2 BCIZ Element Veubeke Element

20

6

Mesh size (N)

Mesh size (N) Fig. 8. Relative error Z in transverse displacement at point C versus mesh size for simply supported square plate and mesh type A.

4

ð27Þ

Is plotted as shown in Figs. 8–11, where wCEX denotes the exact transverse displacement at point C and wC the transverse displacement from the numerical solution. It can be observed that the convergence rate for the transverse displacement at point C obtained using the current element from the new discrete scheme is superior to that of BCIZ and Veubeke elements in all cases.

15 10 5 0 0

2

4

6

8

10

12

14

16

18

Mesh size (N) Fig. 9. Relative error Z in transverse displacement at point C versus mesh size for clamped square plate and mesh type A.

20 18

Element 2 BCIZ Element Veubeke Element

16 Relative error %

14 12 10 8 6 4 2 0 0

2

4

6

8 10 12 Mesh size (N)

14

16

18

Fig. 10. Relative error Z in transverse displacement at point C versus mesh size for Simply supported square plate and mesh type B.

7. Conclusions The weak forms, Eq. (4) or Eq. (7), present a more general form for some discrete schemes, including FEM and the method given in this paper. Using internal moments as intermediate variables, boundary conditions can be introduced easily in weak forms presented with Eq. (12). FEM can be pushed forward from the weak forms, and two reasons (are listed in Section 4) for nonconforming problems in FEM are presented, which may provide an alternative perspective for better understanding of the nonconforming problems in FEM. Based on the given weak forms, a new discrete scheme can be derived in the analysis of the discontinuity of the element interpolation functions. Furthermore, the boundary condition can be easily treated by the use of the equations. Comparing with FEM, the new discretization scheme has clear physical meaning and is not affected by nonconforming problems. The patch test results as shown in Tables 2 and 3 indicate that the current method does not result in loss of accuracy and is not subjected to the irregularity and distortion of finite element mesh. The main source of the error in the current formulation is due to the approximation of the displacement field by the element interpolation function. It should be noted that the new discrete scheme needs further studies to investigate whether the discontinuity of piecewise polynomials across elements causes the loss of accuracy in solution.

Z.X. Tian, D. Wang / Finite Elements in Analysis and Design 47 (2011) 200–207

References [1] J.N. Reddy, in: Energy Principles and Variational Methods in Applied Mechanics, Wiley, New York, 2002. [2] R.J. Melosh, A stiffness matrix for the analysis of thin plates in bending, J. Aeronaut. Sci. 28 (1) (1961) 34–42. [3] B. Fraeijs de Veubeke, A conforming finite element for plate bending, Int. J. Solids Struct. 4 (1968) 95–108. [4] B.M. Irons, A conforming quadratic triangular element for plate bending, Int. J. Numer. Methods Eng. 1 (1) (1969) 29–45. [5] G.P. Bazeley, Y.K. Cheung, B.M. Irons, O.C. Zienkiewicz, Triangular elements in plate bending—conforming and nonconforming solutions, in: Proceedings of the Conference on Matrix Methods in Structural Mechanics, Ohio, October 1965, pp. 547–576. [6] O.C. Zienkiewicz, in: The Finite Element Methods in Engineering Science, McGraw-Hill, New York, 1977. [7] O.C. Zienkiewicz, R.L. Taylor, J.M. Too, Reduced integration technique in general analysis of plates and shells, Int. J. Numer. Methods Eng. 3 (2) (1971) 275–290. [8] E.D.L. Pugh, E. Hinton, O.C. Zienkiewicz, Study of quadrilateral plate bending elements with reduced integration, Int. J. Numer. Methods Eng. 12 (7) (1978) 1059–1079. [9] C. Zienkiewicz, E. Hinton, Reduced integration, function smoothing and nonconformity in finite element analysis, J. Franklin Inst. 302 (5-6) (1976) 443–461. [10] T.J.R. Hughes, R.L. Taylor, W. Kanoknukulcha, A simple and efficient finite element for plate bending, Int. J. Numer. Methods Eng. 11 (10) (1977) 1529–1543. [11] T. Belytschko, C.S. Tsay, Stabilization procedure for quadrilateral plate element with one-point quadrature, Int. J. Numer. Methods Eng. 19 (3) (1983) 405–419. [12] K.J. Bathe, E.N. Dvorkin, A four-node plate bending element based on Mindlin/ Reissner plate theory and a mixed interpolation, Int. J. Numer. Methods Eng. 21 (2) (1985) 367–383. [13] T.H.H. Pian, Derivation of element stiffness matrices by assumed stress distributions, AIAA J. 2 (7) (1964) 1333–1336. [14] B.M. Irons, Engineering application of numerical integration in stiffness method, AIAA J. 14 (1966) 2035–2037.

207

[15] L.M. Tang, W.J. Chen, Y.X. Liu, String net function approximation and quasiconforming technique, in: Proceedings of the International Symposium on Hybrid and Mixed Finite Element Methods, Atlanta, GA, April 2, 1981, pp. 8–10. [16] T.H.H. Pian, D.P. Chen, Alternative ways for formulation of hybrid stress elements, Int. J. Numer. Methods Eng. 18 (1982) 1679–1684. [17] Y.Q. Long, Advances in variational principles in China, in: Proceedings of the 2nd International Conference on Computing in Civil Engineering, China, 1985, pp. 1207–1215. [18] W.J. Chen, Y.K. Cheung, A new approach for the hybrid element method, Int. J. Numer. Methods Eng. 24 (1987) 1697–1709. [19] O.C. Zienkiewicz, T.L. Taylor, J.Z. Zhu, in: The Finite Element Method: Its Basis & Fundamentals, 6th ed., Elsevier Pte Ltd., Singapore, 2005. [20] J.A. Essers, M. Delanaye, P. Rogiest, Upwind-biased finite-volume technique solving Navier–Stokes equations on irregular meshes, AIAA J. 33 (5) (1995) 833–842. [21] C. Debiez, A. Dervieux, K. Mer, B. Nkonga, Computational of unsteady flows with mixed finite volume finite element upwind methods, Int. J. Numer. Methods Fluids 27 (1998) 193–206. [22] P.C. Sames, T.E. Schellin, S. Muzaferija, M. Peric, Application of a two-fluid finite volume method to ship slamming, J. Offshore Mech. Arct. 121 (1) (1999) 44–56. [23] J.S. Wang, H.G. Ni, A high resolution finite volume method for solving shallow water equations, J. Hydrodyn. Ser. B 1 (2000) 35–41. [24] E. Onate, M. Cervera, O.C. Zienkiewicz, Finite volume format for structural mechanics, Int. J. Numer. Methods Eng. 37 (1994) 181–201. [25] M.A. Wheel, A finite volume method for analysing the bending deformation of thick and thin plates, Comput. Methods Appl. Mech. Eng. 147 (1997) 199–208. [26] N. Fallah, A cell vertex and cell centred finite volume method for plate bending analysis, Comput. Methods Appl. Mech. Eng. 193 (2004) 3457–3470. [27] P. Wenke, M.A. Wheel, A finite volume method for solid mechanics incorporating rotational degrees of freedom, Comput. Struct. 81 (2003) 321–329. [28] J. Robinson, Element evaluation. A set of assessment points and standard tests, in: Proceedings of the F.E.M. in the Commercial Environment, vol. 1, October 1978, pp. 217–248.