Two generalized conforming quadrilateral Mindlin–Reissner plate elements based on the displacement function

Finite Elements in Analysis and Design 99 (2015) 24–38

Contents lists available at ScienceDirect

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

Two generalized conforming quadrilateral Mindlin–Reissner plate elements based on the displacement function Yan Shang a,c, Song Cen a,c,d,n, Chen-Feng Li b, Xiang-Rong Fu e a

Department of Engineering Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China Zienkiewicz Centre for Computational Engineering & Energy Safety Research Institute, College of Engineering, Swansea University, Swansea SA2 8PP, UK c High Performance Computing Center, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China d Key Laboratory of Applied Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China e College of Water Conservancy & Civil Engineering, China Agricultural University, Beijing 100083, China b

art ic l e i nf o

a b s t r a c t

Article history: Received 26 June 2014 Received in revised form 7 January 2015 Accepted 16 January 2015 Available online 27 February 2015

This work presents two 4-node, 12-DOF quadrilateral displacement-based finite elements for analysis of the Mindlin–Reissner plate. Derived from the fundamental analytical solutions of the displacement function F, the deflection and rotation fields of the proposed elements satisfy a priori all related governing equations. The unknown coefficients are determined through the generalized conforming element method, a relaxed and rational conforming approach. The resulting elements perform like nonconforming elements on a coarse mesh, and with mesh refinement they converge as conforming elements. Numerical benchmarks demonstrate that the new elements are insensitive to mesh distortion and free of shear locking, and can provide satisfactory results for most cases. & 2015 Elsevier B.V. All rights reserved.

Keywords: Finite element Mindlin–Reissner plate Displacement function Generalized conforming Mesh distortion

1. Introduction As a major class of finite element models, Mindlin–Reissner plate bending elements have been extensively studied over the past few decades and numerous successful models were proposed [1–15]. Playing a central role in practical engineering analysis, the Mindlin– Reissner plate element is continuously being improved for better performance and to meet specific requirement. Among others, some of the latest developments include the NIPE formulations derived from the assumed-strain technique [16], the models with continuous displacement and discontinuous rotations [17,18], the element based on the high-order linked interpolation [19], the alternative alpha finite element method with discrete shear gap technique [20], and so on. However, two major challenges remain outstanding: (1) the performance of the element cannot be guaranteed when the mesh is severely distorted; and (2) the precisions for stress/resultant results are relatively lower than those for displacements. A series of new developments have been reported recently, which make it possible to address the aforementioned challenges. Liu et al. [21] proposed the smoothed FEM (SFEM), which integrates the strain smoothing technique into the conventional FEM and formulates different smoothed FEM models, including the cell-based SFEM

n Corresponding author at: Department of Engineering Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China. E-mail address: [email protected] (S. Cen).

http://dx.doi.org/10.1016/j.finel.2015.01.012 0168-874X/& 2015 Elsevier B.V. All rights reserved.

[22,23], the edge-based SFEM [24] and the node-based SFEM [25]. The plate elements based on these SFEM models [26,27] are insensitive to mesh distortion and free of shear locking. Ribatić et al. [28] constructed two 9-node distortion-immune displacement-based quadrilateral thick plate elements by using bubble parameters. Cen et al. [29] proposed a hybrid displacement function (HDF) method for formulating Mindlin–Reissner plate elements, where they took the displacement function F [30] to derive the stress trial functions that satisfy all governing equations and employed the Timoshenko’s beam formulae to determine boundary displacement modes. The HDF element enables superior precision for both displacements and resultants, even in cases where a severely distorted mesh containing concave quadrilateral or degenerated triangular elements is used. By combining the displacement function F [29,30] with the genera lized conforming element method [1], this work develops two highperformance displacement-based Mindlin–Reissner plate elements. The concept of the generalized conforming element was first proposed by Long et al. [31], and since then has been successfully applied to plane problems [32,33], plate problems [34–37] and shell problems [38] as well. In this concept, the displacement fields are usually determined by the relaxed compatibility requirements, i.e., the generalized conforming conditions. The new generalized conforming elements are derived via three steps. First, the deflection and rotation fields within the element are assumed according to the fundamental analytical solutions of the displacement function F, and the corresponding unknown coefficients

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

are introduced. Next, the relations between these unknown coefficients and the element’s nodal displacement DOFs are determined by employing a set of generalized compatibility conditions, from which the deflection and rotation fields can be expressed as functions of the element’s nodal displacement DOFs. Finally, the new elements can be constructed following the principle of minimum potential energy. These two quadrilateral elements are denoted as GCP4-12α and GCP4-14α, respectively, indicating that they are “Generalized Conforming Plate element with 4 nodes and 12 or 14 displacement coefficients”. Since the displacement components derived from the displacement function F can satisfy all related governing equations, these new elements also possess the advantages from analytical methods, as the HDF elements [29] do. The performance of the new elements is assessed with standard benchmark examples. Numerical results show that these two new elements are free of shear locking, insensitive to mesh distortions, and can provide satisfactory results for most cases. Only a small deviation exists for thick plates under soft simply supported boundary conditions due to interelement incompatibility.

2. The formulations of new elements 2.1. The displacement function F For a Mindlin–Reissner plate (see Fig. 1), the transverse deflection w and the rotations ψx, ψy can be expressed by the displacement ψx

z,w

ψy

∂ w/∂ x x,u

D w ¼ F  ∇2 F; C

ψx ¼

y,v

Mxy

Mx

ð1Þ

ð2Þ

Substitution of Eq. (1) into the strain–displacement relations yields κxx ¼ 

∂2 F ; ∂x2

γ xz ¼ 

D ∂
2  ∇ F ; C ∂x

κ yy ¼ 

∂2 F ; ∂y2

γ yz ¼ 

∂2 F ; κxy ¼  2 ∂x∂y

ð3Þ

D ∂
2  ∇ F ; C ∂y

ð4Þ

where D and C denote the bending and shear stiffness of the plate. The solutions of the deflection, rotations and strains can be readily derived from the fundamental analytical solutions of displacement function F [29]. For the general (homogeneous) part of F, the displacement components and strains are listed in Table 1. For the particular part corresponding to a plate subjected to a uniformly distributed transverse load q, the displacement components are  q
2 9 8 n9 8 q
4 x þ y4  4C x þ y2 > w > > > > = < 48D < = > n q 3 ; ð5Þ un ¼ ψ x ¼ 12Dx > > > > : ψn > ; > q 3 ; : y 12Dy

A3

7

B3

3

A4

B2

My

8

6

Tx

Gauss points

A2

1

My Ty Mid-surface (xoy plane)

Fig. 1. The definitions of the displacements and resultants for a Mindlin–Reissner plate.

nodes

Mid-side points

B4

Mxy

h Mxy

∂F ; ∂y

D∇2 ∇2 F ¼ q:

Mx

x

ψy ¼

Mxy

Ty Tx

∂F ; ∂x

in which F should satisfy the equation:

z,w

y

z

function F [30]:

4

∂w/∂ y

25

A1

5

B1

2

Fig. 2. The generalizing conforming quadrilateral plate element with 4 nodes.

Table 1 Fourteen fundamental analytical solutions for the general part F0 of the displacement function and the corresponding deflections, rotations and strains. i

1

2

3

4

F 0i

1

x

y

x2

xy

y2

x3

x2y

xy2

w0i

1

x

y

x2  2D/C

xy

y2  2D/C

x3  6xD/C

x2y  2yD/C

xy2  2xD/C

ψ 0xi

0

1

0

2x

y

0

3x2

2xy

y2

ψ 0yi

0

0

1

0

x

2y

0

x2

2xy

κ 0xxi

0

0

0

2

0

0

 6x

 2y

0

κ 0yyi

0

0

0

0

0

2

0

0

 2x

κ 0xyi

0

0

0

0

2

0

0

 4x

 4y

γ 0xzi

0

0

0

0

0

0

 6D/C

0

 2D/C

γ 0yzi

0

0

0

0

0

0

0

 2D/C

0

i

10 y3

11 x3y

12 xy3

13 x4  y4

14 6x2y2  x4  y4

w0i

y3  6yD/C

x3y  6xyD/C

xy3  6xyD/C

x4  y4  (12x2  12y2)D/C

6x2y2  x4  y4

ψ 0xi

0

3x2y

y3

4x3

12xy2  4x3

3y2

x3

3xy2

 4y3

F 0i

ψ 0yi κ 0xxi κ 0yyi κ 0xyi γ 0xzi γ 0yzi

0

 6xy

0

5

12x2y 4y3 2

 12x 2

 12y2 þ 12x2  12x2 þ12y2

 6y

0

 6xy

12y

0

 6x2

 6y2

0

 48xy

0

 6yD/C

 6yD/C

 24xD/C

0

 6D/C

 6xD/C

 6xD/C

24yD/C

0

6

7

8

9

26

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

4-node quadrilateral element, as shown in Fig. 2, the displacement function F is set as the linear combination of its fundamental analytical solutions:

and the strains are 8 n 9 8 q 29 κxx > > > >  4Dx > > > > > > > > > > > > q 2> > > > > κnyy >  y > > > = = < 4D > < n n 0 : E ¼ κ xy ¼ > > > > > > > > >  qx > > γ nxz > > > > > 2C > > > > > > ; > > : γn > ; :  qy > yz

ð6Þ

F ¼ F0 þ Fn ¼

k X

F 0i αi þF n ;

ð7Þ

i¼1

2C

2.2. The procedure for constructing new generalized conforming elements Unlike the traditional generalized conforming elements or nonconforming elements, the deflection and rotation fields of the new elements are assumed according to the fundamental analytical solutions of the displacement function F. For a generalized conforming

where k is the number of the fundamental analytical solutions for the general part of F0, and αi (i¼1 k) are the corresponding unknown coefficients. Substituting Eq. (7) into Eq. (1), the deflection and rotation fields within the element are obtained as

w¼

k X

w0i αi þ wn ;

i¼1

Fig. 3. The model and mesh for the patch test.

ψx ¼

k X i¼1

ψ 0xi αi þ ψ nx ;

ψy ¼

k X i¼1

ψ 0yi αi þ ψ ny :

ð8Þ

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

The potential energy functional of a Mindlin–Reissner plate element can be expressed in the following matrix form [1]:

The displacements can be expressed in the following matrix form: 2 u¼

w01 6 0 6 ψ x1 4 ψ 0y1

w02 ψ 0x2 ψ 0y2

w03 ψ 0x3 ψ 0y3

⋯ ⋯ ⋯

8 9 α > > 1> 3> > 8 n9 0 > wk > w > > > α2 > > > = > = < < 0 7 ψ xk 7 α3 þ ψ nx ¼ Uα þ un : 5> > > > : ψn > ; > ψ 0yk > > > ⋮ > y > > > ; :α >

1 T Π eP ¼ ∬Ae ET AEdxdy ∬Ae f udxdyþ 2 ð9Þ

e

λα þ U ¼ Γq ; where h qe ¼ w 1

ψ x1

ð10Þ

Mesh N  N

w2

ψ x2

ψ y2

w3

ψ x3

ψ y3

w4

ψ x4

ψ y4

iT

22

ð19Þ

R dds;

88

16  16

32  32

Reference

/100D) 0.1259 0.1259 0.1251 0.1238 0.1263 0.1263 0.1264

0.1264 0.1264 0.1262 0.1259 0.1265 0.1265 0.1265

0.1265 0.1265 0.1264 0.1265 0.1265 0.1265 0.1265

0.1265 0.1265 0.1265 0.1265 0.1265 0.1265 0.1265

0.1265

Central moment Mc/(qL2/10) GCP4-12α 0.2234 0.2283 GCP4-14α 0.2234 0.2283 MITC4 [2,27] 0.1890 0.2196 Q4BL [10] 0.2179 0.2261 ARS-Q12 [14] 0.2869 0.2433 AC-MQ4 [37] 0.2712 0.2407 MISC2 [27] 0.1976 0.2218

0.2290 0.2290 0.2267 0.2283 0.2326 0.2321 0.2273

0.2290 0.2290 0.2285 0.2289 0.2299 0.2298 0.2286

0.2291 0.2291 0.2289 0.2290 0.2293 0.2292 0.2289

0.2291

Central deflection wc/(qL GCP4-12α 0.1227 GCP4-14α 0.1227 MITC4 [2,27] 0.1211 Q4BL [10] 0.1116 ARS-Q12 [14] 0.1245 AC-MQ4 [37] 0.1245 MISC2 [27] 0.1266

ð11Þ

ψ y1

T

Seσ

Table 2 The dimensionless results of central deflection wc/(qL4/100D) and moment Mc/(qL2/ 10) of the clamped square plate, thin plate case (h/L ¼0.001).

The detailed components of matrices U and E in the above equations are listed in Table 1. In order to solve these unknown coefficients αi (i ¼1  k), it is necessary to introduce k generalized conforming conditions. After employing an appropriate set of generalized conforming conditions, the relations between the coefficients αi (i ¼1  k) and the element’s nodal displacement DOFs are obtained as n

Z

where A is the elasticity matrix as defined in Eq. (20), f is the distributed transverse load vector within the plate as defined in Eq. (21), R is the boundary load vector at the element’s edge as

k

Then, the corresponding strain field can be obtained as 9 2 38 9 8 κ xx1 κ xx2 κxx3 ⋯ κxxk > α1 > > κnxx > > > > > > > > n > > > > 6 κyy1 κyy2 κyy3 ⋯ κ 7> yyk 7> > > > > α2 > > κ yy > 6 6 7< = < n = 6 7 κ ¼ Eα þ En : E ¼ 6 κ xy1 κ xy2 κ xy3 ⋯ κ xyk 7 α3 þ xy > > > n > > > > > 6γ 7 > > > > γ ⋮ > > xz > 4 xz1 γ xz2 γ xz3 ⋯ γ xzk 5> > > > > > > > > γ yz1 γ yz2 γ yz3 ⋯ γ yzk : αk ; : γ nyz ;

27

:

44 4

ð12Þ n

The matrices λ, U and Γ have different expressions for different sets of generalized conforming conditions. The details for this matter will be discussed in Section 2.3. Following Eq. (11), α can be expressed in terms of qe: α ¼ Lλ qe  Lu ;

ð13Þ

where Lλ ¼ λ  1 Γ;

L u ¼ λ  1 Un :

ð14Þ

Then, substitution of Eq. (13) into Eqs. (9) and (10) yields
 ð15Þ u ¼ U Lλ qe  Lu þun ¼ Nqe þ Nn ;
 E ¼ E Lλ qe  Lu þ En ¼ Bqe þ Bn ;

ð16Þ

in which N ¼ ULλ ;

Nn ¼ un  ULu ;

ð17Þ

B ¼ ELλ ;

Bn ¼ En  ELu :

ð18Þ

Table 3 The dimensionless results of central deflection wc/(qL4/100D) and moment Mc/(qL2/ 10) of the clamped square plate, thick plate case (h/L ¼0.1). Mesh N  N

22

88

16  16

32  32

Reference

/100D) 0.1499 0.1499 0.1488 0.1491 0.1548 0.1477 0.1500

0.1501 0.1501 0.1500 0.1501 0.1515 0.1494 0.1503

0.1502 0.1502 0.1504 0.1504 0.1507 0.1502 0.1504

0.1502 0.1502 0.1504 0.1504 0.1505 0.1504 0.1505

0.1499

Central moment Mc/(qL2/10) GCP4-12α 0.2350 0.2306 GCP4-14α 0.2350 0.2306 MITC4 [2,27] 0.1898 0.2219 Q4BL [10] 0.2186 0.2288 ARS-Q12 [14] 0.3061 0.2527 AC-MQ4 [37] 0.2813 0.2494 MISC2 [27] 0.1982 0.2241

0.2310 0.2310 0.2295 0.2313 0.2376 0.2373 0.2300

0.2315 0.2315 0.2314 0.2318 0.2334 0.2334 0.2315

0.2317 0.2317 0.2318 0.2320 0.2324 0.2324 0.2319

0.231

Central deflection wc/(qL GCP4-12α 0.1494 GCP4-14α 0.1494 MITC4 [2,27] 0.1431 Q4BL [10] 0.1427 ARS-Q12 [14] 0.1678 AC-MQ4 [37] 0.1474 MISC2 [27] 0.1483

44 4

L /3

y

L /3

ψx = 0

L /3

x c

ψy = 0

L /3

Fig. 4. A quarter of square plate (c is the central point of plate). (a) 4×4 regular mesh and (b) 4×4 distorted mesh.

28

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

defined in Eq. (22), and d is the boundary displacement along the boundary Seσ as defined in Eq. (23). 2

D 6 μD 6 6 A¼6 6 0 6 4 0 0

μD

0

0

D

0

0

0 0

1μ 2 D

0

0 C

0

0

0

0

3

Mesh N  N

07 7 7 07 7; 7 05 C

ð20Þ

Table 4 The dimensionless results of central deflection wc/(qL4/100D) and moment Mc/(qL2/ 10) of the SS1 square plate, thin plate case (h/L¼ 0.001). Mesh N  N

22

44

Central deflection wc/(qL4/100D) GCP4-12α 0.4063 0.4062 GCP4-14α 0.4063 0.4062 Q4BL [10] 0.4593 0.4292 ARS-Q12 [14] 0.4045 0.4060 AC-MQ4 [37] 0.4033 0.4061 Central moment Mc/(qL2/10) GCP4-12α 0.4774 0.4788 GCP4-14α 0.4774 0.4788 Q4BL [10] 0.5649 0.5010 ARS-Q12 [14] 0.5005 0.4839 AC-MQ4 [37] 0.5897 0.5050

88

16  16

32  32

Reference

0.4062 0.4062 0.4164 0.4062 0.4062

0.4062 0.4062 0.4110 0.4062 0.4062

0.4062 0.4062 0.4086 0.4062 0.4062

0.4062

0.4789 0.4789 0.4876 0.4800 0.4850

0.4789 0.4789 0.4830 0.4792 0.4799

0.4789 0.4789 0.4809 0.4789 0.4790

22

0.4789

88

16  16

32  32

Reference

/100D) 0.4500 0.4500 0.4727 0.4419 0.4437

0.4532 0.4532 0.4644 0.4544 0.4543

0.4542 0.4542 0.4624 0.4596 0.4595

0.4544 0.4544 0.4618 0.4612 0.4611

0.4617

Central moment Mc/(qL2/10) GCP4-12α 0.4895 0.4976 GCP4-14α 0.4895 0.4976 Q4BL [10] 0.5694 0.5169 ARS-Q12 [14] 0.5206 0.5087 AC-MQ4 [37] 0.6115 0.5236

0.5010 0.5010 0.5112 0.5081 0.5122

0.5025 0.5025 0.5100 0.5091 0.5101

0.5029 0.5029 0.5100 0.5094 0.5097

0.5096

Central deflection wc/(qL GCP4-12α 0.4430 GCP4-14α 0.4430 Q4BL [10] 0.4957 ARS-Q12 [14] 0.4280 AC-MQ4 [37] 0.4358

44

88

16  16

32  32

Reference

Central deflection wc/(qL4/100D) GCP4-12α 0.4225 0.4256 GCP4-14α 0.4225 0.4256 MITC4 [2,27] 0.4190 0.4255 Q4BL [10] 0.4269 0.4274 ARS-Q12 [14] 0.4228 0.4255 AC-MQ4 [37] 0.4139 0.4206 MISC2 [27] 0.4285 0.4277

0.4268 0.4268 0.4268 0.4273 0.4267 0.4251 0.4274

0.4272 0.4272 0.4272 0.4273 0.4271 0.4267 0.4273

0.4273 0.4273 0.4273 0.4273 0.4272 0.4271 0.4273

0.4273

Central moment Mc/(qL2/10) GCP4-12α 0.4810 0.4767 GCP4-14α 0.4810 0.4767 MITC4 [2,27] 0.4075 0.4612 Q4BL [10] 0.4716 0.4773 ARS-Q12 [14] 0.5223 0.4941 AC-MQ4 [37] 0.5419 0.5028 MISC2 [27] 0.4172 0.4637

0.4778 0.4778 0.4745 0.4785 0.4834 0.4862 0.4751

0.4785 0.4785 0.4778 0.4788 0.4801 0.4808 0.4779

0.4788 0.4788 0.4786 0.4788 0.4792 0.4794 0.4786

0.4789

 f¼ q

Table 5 The dimensionless results of central deflection wc/(qL4/100D) and moment Mc/(qL2/ 10) of the SS1 square plate, thick plate case (h/L ¼ 0.1). Mesh N  N

Table 7 The dimensionless results of central deflection wc/(qL4/100D) and moment Mc/(qL2/ 10) of the SS2 square plate, thick plate case (h/L ¼0.1).

4

0

h R ¼ T

22

0

T

Mn

44

;

ð21Þ M ns

8 9 2  1 > =  6 d ¼ ψn  ¼ 4 > :ψ > ; s e Sσ

iT

;

l m

ð22Þ 38 9  w  > < > =  7 m 5 ψ x  ¼ Luj Seσ ; > :ψ > ; l y e

ð23Þ

Sσ

where l, m denote the direction cosines of the element boundaries’ outer normal n. By substituting Eqs. (15) and (16) into Eq. (19) and applying the principle of minimum potential energy, ∂Π eP ¼ ∬Ae BT ABdxdyqe þ ∬Ae BT ABn dxdy  ∬Ae NT fdxdy ∂qe Z NjTSeσ LT Rds ¼ 0; þ Sσ e

ð24Þ

the following relation can be obtained: Kqe ¼ pe ;

ð25Þ

where K is the element stiffness matrix K ¼ ∬Ae BT ABdxdy;

ð26Þ

e

Table 6 The dimensionless results of central deflection wc/(qL4/100D) and moment Mc/(qL2/ 10) of the SS2 square plate, thin plate case (h/L¼ 0.001). Mesh N  N

22

88

16  16

32  32

Reference

/100D) 0.4062 0.4062 0.4041 0.4058 0.4060 0.4062 0.4077

0.4062 0.4062 0.4057 0.4062 0.4062 0.4062 0.4066

0.4062 0.4062 0.4061 0.4062 0.4062 0.4062 0.4063

0.4062 0.4062 0.4062 0.4062 0.4062 0.4062 0.4063

0.4062

Central moment Mc/(qL2/10) GCP4-12α 0.4786 0.4788 GCP4-14α 0.4786 0.4788 MITC4 [2,27] 0.4075 0.4612 Q4BL [10] 0.4712 0.4773 ARS-Q12 [14] 0.5009 0.4839 AC-MQ4 [37] 0.5106 0.4872 MISC2 [27] 0.4171 0.4637

0.4789 0.4789 0.4745 0.4784 0.4801 0.4810 0.4751

0.4789 0.4789 0.4778 0.4788 0.4792 0.4794 0.4779

0.4789 0.4789 0.4786 0.4788 0.4789 0.4790 0.4786

0.4789

Central deflection wc/(qL GCP4-12α 0.4051 GCP4-14α 0.4051 MITC4 [2,27] 0.3969 Q4BL [10] 0.4305 ARS-Q12 [14] 0.4045 AC-MQ4 [37] 0.4052 MISC2 [27] 0.4123

44 4

and p is the element nodal equivalent load vector Z pe ¼ ∬Ae NT fdxdy  ∬Ae BT ABn dxdy NjTSeσ LT Rds; Seσ

ð27Þ

where Nj Seσ is the value of N along the boundary Seσ . The unknown coefficient α can be derived from Eq. (13) after solving for qe. Then the displacements and the strains at any point of the element can be determined by substituting its Cartesian coordinates into Eqs. (9) and (10). The above derivation shows the general formulations for constructing new generalized conforming elements. If different generalized conforming conditions are employed for deriving the elements, the detailed expressions of related matrices (see from Eqs. (14)–(18)) will also be different. 2.3. The formulations of the new generalized conforming plate (GCP) elements In this section, two new generalized conforming plate (GCP) elements are constructed following the procedure described in Section 2.2. According to the number of the fundamental analytical

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

solutions of the displacement function F being used, these two elements are denoted by GCP4-12α and GCP4-14α, respectively. 2.3.1. The element GCP4-12α As shown in Fig. 2, for the element GCP4-12α, the displacement function F (see Eq. (7)) is assumed to be the linear combination of its first 12 terms of fundamental analytical solutions (see Table 1). Then, k in Eqs. (7–10) equals 12. To determine these coefficients αi (i¼1–12), 12 generalized conforming conditions must be employed. First, the following nodal compatibility conditions for deflection are considered:   w xj ; yj ¼ wj ; ðj ¼ 1 4Þ; ð28Þ where (xj, yj), (j¼1–4), are the Cartesian coordinates of the node j, w(xj, yj) is calculated by using Eq. (9), and as one of the element’s DOFs, wj is the deflection at the node j. Next, the point compatibility conditions are considered for normal rotations at two Gauss pffiffiffi points Ai and Bi (the local parametric coordinates are 7 1= 3, respectively) along each edge:
 ð29Þ ψ n xk ; yk ¼ ψ nk ; ðk ¼ A1 ; B1 ; A2 ; B2 ; A3 ; B3 ; A4 ; B4 Þ; where ψn(xk, yk) are also calculated by substituting the Cartesian coordinates into Eq. (9), and ψ nk are determined by the boundary normal rotations ψ n varying linearly along the edge ij: ψ nij ¼ ð1  sÞψ ni þ sψ nj ;

ðs ¼ 0  1Þ:

29

where ψ ni and ψ nj are the normal rotations at the node i, j along the edge ij. Substituting Eqs. (7–10) into Eqs. (28–30), the matrices in Eq. (11) can be obtained for the element GCP4-12α. Specifically, the matrix λ can be written as: " λ¼

W

#

φn

;

φn ¼ Tn Φ;

ð31Þ

in which 2


 w01 x1 ; y1  6 0
6 w1 x2 ; y2 6  W ¼ 6 0
6 w1 x3 ; y3 4
 w01 x4 ; y4


 w02 x1 ; y1
 w02 x2 ; y2
 w02 x3 ; y3
 w02 x4 ; y4


 w03 x1 ; y1
 w03 x2 ; y2
 w03 x3 ; y3
 w03 x4 ; y4

⋯ ⋯ ⋯ ⋯


3 w012 x1 ; y1
7 w012 x2 ; y2 7
7 7; w012 x3 ; y3 7 5
 w012 x4 ; y4 ð32Þ


 ψ 0x1 xA1 ; yA1  6 0
6 ψ y1 xA1 ; yA1 6 6 ⋮ Φ¼6  6 0
6 ψ x1 xB4 ; yB4 4
 ψ 0y1 xB4 ; yB4 2


 ψ 0x2 xA1 ; yA1
 ψ 0y2 xA1 ; yA1 ⋮


 ψ 0x2 xB4 ; yB4
 ψ 0y2 xB4 ; yB4

⋯ ⋯ ⋱ ⋯ ⋯


3 ψ 0x12 xA1 ; yA1
7 ψ 0y12 xA1 ; yA1 7 7 7 ⋮ ;
7 7 ψ 0x12 xB4 ; yB4 7 5
 ψ 0y12 xB4 ; yB4

ð33Þ

ð30Þ

Fig. 5. Errors of the central deflections and moments for clamped square plates in regular mesh. (a) h/L=0.001 (thin plate case) and (b) h/L=0.1 (thick plate case).

30

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

2 6 6 Tn ¼ 6 4

2

3

T1

"

7 7 7; 5

T2 T3

Ti ¼

lij


3 ψ nx xA1 ; yA1
7 6 n 6 ψ y xA1 ; yA1 7 6 7 6 7 Φn ¼ 6
⋮ 7:  6 n 7 6 ψ x xB4 ; yB4 7 4
5 n ψ y xB4 ; yB4

#

mij lij

mij

;

ð34Þ

T4



!



! where i ¼ 1; 2; 3; 4 ; j ¼ 2; 3; 4; 1 , and lij and mij denote the direction cosines of the outer normal of the element edge ij. The matrix Un can be written as " n# W Un ¼ ð35Þ ; φnn ¼ Tn Φn ; φnn
3 wn x1 ; y1
7 6 n 6 w x2 ; y2 7
7 Wn ¼ 6 6 wn x3 ; y3 7; 4
5 wn x4 ; y4

The matrix Γ is given by " # Γw Γ¼ ; Γφ

2

2

0 6 60 6 6 60 6 6 60 6 Γφ ¼ 6 60 6 6 60 6 6 60 6 4 0

1 ð1  sA1 Þy12  l12

2 6 6 Γw ¼ 6 4

ð36Þ

ð38Þ

3

I13

7 7 7; 5

I13 I13

 I13 ¼ 1

0

 0 ;

ð39Þ

I13

0

1  l12 sA1 y12

 l12 ð1  sB1 Þy12

1 ð1 sA1 Þx12 l12 1 ð1  sB1 Þx12 l12

0

 l12 sB1 y12

0

0

0

1 ð1  sA2 Þy23  l23

1

ð37Þ

1

1 s x l12 A1 12 1 s x l12 B1 12

0

0

0

0

0

0

0

0

0

0

0

1  l23 sA2 y23

0

0

0

 l23 sB2 y23

1 s x l23 A2 23 1 s x l23 B2 23

0

0

0

1  l34 sA3 y34

0

0

0

 l23 ð1 sB2 Þy23

1 ð1  sA2 Þx23 l23 1 ð1  sB2 Þx23 l23

0

0

0

0

0

0

1  l34 ð1 sA3 Þy34

1

1

0

0

0

0

0

1  l34 ð1  sB3 Þy34

1 ð1  sA3 Þx34 l34 1 ð1 sB3 Þx34 l34

0

1  l34 sB3 y34

1  l41 sA4 y41

1

s x l41 A4 41

0

0

0

0

0

0

0

1  l41 ð1  sA4 Þy41

1  l41 sB4 y41

1 s x l41 B4 41

0

0

0

0

0

0

0

1  l41 ð1  sA4 Þy41

0

0

3

7 7 7 7 7 0 7 7 7 0 7 7; 1 7 s x 7 l34 A3 34 7 1 7 s x B3 34 7 l34 7 1 7 ð 1 s Þx 41 A4 7 l41 5 1 ð1 sA4 Þx41 l41 0

Fig. 6. Errors of the central deflections and moments for SS1 square plates in regular mesh. (a) h/L=0.001 (thin plate case) and (b) h/L=0.1 (thick plate case).

ð40Þ

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

where xij ¼ xi  xj ; yij ¼ yi  yj , and lij is the length of edge ij. Then, by substituting the matrices λ, Un and Γ (see Eqs. (31–40)) into the equations discussed in Section 2.2, the expressions of the stiffness matrix and the nodal equivalent load vector for the element GCP4-12α can be obtained.

31

beam [29]: wij ¼ ½1  s þð1 2δij ÞZ 3 wi þ ½s  ð1  2δij ÞZ 3 wj lij lij þ ½Z 2 þ ð1  2δij ÞZ 3 ψ si  ½Z 2  ð1  2δij ÞZ 3 ψ sj ; 2 2

ð42Þ

in which 2.3.2. The element GCP4-14α For the element GCP4-14α (see Fig. 2), the displacement function F is formed by the first 14 terms of its fundamental analytical solutions (see Table 1). To determine the unknown coefficients αi (i¼1–14), fourteen generalized conforming conditions are needed. Besides those conditions used for the element GCP4-12α, two additional compatibility conditions are considered: wðx5 ; y5 Þ þwðx7 ; y7 Þ ¼ w5 þ w7 wðx6 ; y6 Þ þwðx8 ; y8 Þ ¼ w6 þ w8 ;

ð41Þ

where (xj, yj) (j¼ 5–8) are the Cartesian coordinates of the mid-side points 5–8 at each element edge, and wj (j¼ 58) are the deflections at these mid-side points. The deflections wj (j¼5–8) are determined by using the locking-free formulae of Timoshenko’s " Wmid ¼



 w01 x5 ; y5 þ w01 x7 ; y7

 w01 x6 ; y6 þ w01 x8 ; y8



 w02 x5 ; y5 þw02 x7 ; y7

 w02 x6 ; y6 þw02 x8 ; y8

Z 2 ¼ sð1  sÞ;

Z 3 ¼ sð1  sÞð1  2sÞ;

δij ¼

6λij ; 1 þ 12λij

λij ¼

D 2

Clij

;

ð43Þ

where D and C denote the bending and shear stiffness, respecti vely; and lij is the edge length. For these mid-side points 5–8, s equals 1/2. Then, by substituting Eqs. (7–10) into these fourteen generalized conforming conditions (see Eqs. (28)–(30) and (41)–(43)), the details of matrices in Eq. (11) can be obtained. Specifically, the matrix λ is given by 2 3 W 6 φ 7 λ ¼ 4 n 5; ð44Þ Wmid where W and φn share the same for as in Eq. (31), and Wmid is obtained as ⋯ ⋯



# w014 x5 ; y5 þ w014 x7 ; y7

 : w014 x6 ; y6 þ w014 x8 ; y8

Fig. 7. Errors of the central deflections and moments for SS2 square plates in regular mesh. (a) h/L=0.001 (thin plate case) and (b) h/L=0.1 (thick plate case).

ð45Þ

32

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

The matrix Un is 2

n

3

W 6 7 Un ¼ 4 φnn 5; n Wmid

ð46Þ

where Wn and φnn are the same as Eq. (35), and Wnmid is obtained as: " n

# w x5 ; y5 þ wn x7 ; y7 n

 : Wmid ¼ ð47Þ wn x6 ; y6 þ wn x8 ; y8

The matrix Γ is 2 3 Γw 6 7 Γ ¼ 4 Γφ 5; Γnw

ð48Þ

where Γw and Γφ are the same as Eq. (38), and Γnw is obtained as "1

Γnw ¼

2 1 2

 18x12

 18y12

1 8x41

1 8y41

1 2 1 2

1 8x12  18x23

1 8y12  18y23

1 2 1 2

 18x34

 18y34

1 8x23

1 8y23

1 2 1 2

1 8x34  18x41

1 8y34  18y41

#

ð49Þ where xij ¼ xi xj ; yij ¼ yi  yj .

Fig. 8. Errors of the central deflections and moments for square plates in distorted mesh. (a) Clamped case; (b) SS1 case and (c) SS2 case.

;

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

33

For the constant shear deformation state, the exact displacement fields are [46]

Thus, the stiffness matrix and the nodal equivalent load vector for the element GCP4-14α can be readily obtained. Since the element GCP4-14α uses 14 terms of fundamental analytical solutions, which guarantees the completeness of the displacement function up to the fourth order, it performs better than the element GCP4-12α in distorted meshes.

w¼

ðx þ yÞ ; 2

ψx ¼

1 ; 2

ψy ¼

1 : 2

ð52Þ

This test is significant only for the thick plate with a very large thickness in which only shear energy works. The element

3. Numerical tests Several benchmark examples are employed to assess the new elements GCP4-12α and GCP4-14α. In order to clearly illustrate the performance, the results from some well-known Mindlin–Reissner plate elements are also provided for comparison, including MITC4 [2,27], Q4BL [10], ARS-Q12 [14], RDKMQ [15], MISC2 [27], AC-MQ4 [37], S4R [39], CRB1 [40], CRB2 [40], S1 [40], DKQ [41], 9βQ4 [42], MiSP4 [43], MMiSP4 [43], MiSP4 þ [44] and DKMQ [45]. 3.1. The patch test The model and the mesh used for the patch test are shown in Fig. 3a together with the parameter settings. Three cases with different span–thickness ratios (2a/h¼1000, 100, 20) are considered, where 2a is the length of the longer edge and h is the plate’s thickness. (a) Displacement loading case In this case, the deflections and rotations at the outer nodes 1–4 are applied to the patch as the boundary conditions, and their values at the inner nodes 5–8 will be evaluated. For the constant bending state, the exact displacement fields are given by [28]
 1 þ x þ 2y þ x2 þ y2 ð1 þ 2xÞ ð2 þ 2yÞ w¼ ; ψy ¼ : ð50Þ ; ψx ¼ 2 2 2 Both GCP4-12α and GCP4-14α can provide exact solutions in all span–thickness ratio cases (2a/h ¼1000, 100, 20). For the constant twisting state, the exact displacement fields are given by [28] ð1 þ x þ2y þ xyÞ w¼ ; 2

ð1 þ yÞ ψx ¼ ; 2

ð2 þ x Þ ψy ¼ : 2

Fig. 10. Results for the sensitivity test to mesh distortion: top: symmetric case; bottom: asymmetric case.

ð51Þ

L

In case 2a/h¼ 1000, both new elements can pass the test exactly. In case 2a/h ¼100, the maximum relative percentage errors of elements GCP4-12α and GCP4-14α are 0.078% and 0.074%, respectively, and increase to 1.73% and 2.05% in case 2a/h ¼20. These errors are mainly caused by the interelement incompatibilities in a coarse mesh. With the refinement of the mesh, these errors will vanish rapidly.

L/2=5

P

Δ

P L

C

Fig. 11. The clamped cantilever plate for sensitivity test to mesh distortion.

L/2=5

C

L/2=5

L/2=5

Δ

L

Δ Fig. 9. The clamped thin square plate for sensitivity test to mesh distortion.

C

Δ

34

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

GCP4-14α can exactly pass this test, while the maximum relative percentage error of GCP4-12α is 6.0%. (b) Constant bending moment loading case (Mn ¼1) Fig. 3b shows this loading case in which the plate is subjected to moments along its all edges. The equivalent nodal forces are also given in Fig. 3b. In this case, the bending moments Mx ¼1 and My ¼1, while other resultants are equal to zero. For all different span–thickness ratios (2a/h¼1000, 100, 20), both GCP4-12α and GCP4-14α can exactly pass this patch test. (c) Constant twisting moment loading case (Mxy ¼1) As shown in Fig. 3c, the plate is subjected to twisting moment along its all edges. The equivalent nodal forces, derived from the assumption that the plate’s edge behaves like a lockingfree Timoshenko’s beam [29], are also given in Fig. 3c. In this case, the twisting moments Mxy ¼ 1, while other resultants are equal to zero. When 2a/h¼ 1000, both two new elements can pass the test exactly. For the ratio 2a/h¼100, the maximum relative percentage errors of the elements GCP4-12α and GCP4-14α are 2.12% and 1.70%, respectively. And for the ratio 2a/h¼ 20, they increase to 3.24% and 2.80%. (d) Constant shear resultant loading case (Tx ¼1) As stated in Section 3.1(a), this test is significant only for the thick plate limit [45]. Fig. 3d shows the cantilever plate, with a very large value of thickness h, is subjected to a distributed line shear force at the free edge (replaced by two concentrated nodal forces). The value of shear force Tx should equal 1. The element GCP4-14α can exactly pass this test, while the maximum relative percentage error of GCP4-12α is 2.62%.

Fig. 12. Results for deviation of the transverse tip displacement wc.

L

E y

D

E=10.92 μ =0.3 h=0.1 L=100

Free L

C Free 60°

A

x

B Fig. 13. The Razzaque’s 601 skew plate.

Fig. 14. Errors of the central deflection wc and moment My of Razzaque’s skew plate.

Radius R=5 E=10.92; μ =0.3; Uniform loading: q=1 Displacement BCs: w =0 (SS1) along AB w =ψn= ψs=0 (Clamped) along AB

y

y A

A

ψx=0

ψx=0

x

x C

ψy=0

B

C

ψy=0

B

Fig. 15. A quarter of circular plate.

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

35

be eliminated by using a refiner mesh. Since the element GCP414α employs two more trial functions, it indeed exhibits better performance. As described in previous sections, these two elements belong to the nonconforming model in a coarse mesh, and will converge into conforming models when the mesh is refined. Thus, the convergence can still be guaranteed.

(e) Some discussions From above tests, it can be observed that, for the thin plate (2a/ h¼1000), both new elements can provide exact solutions. However, small errors may appear when the thickness increases for constant twisting moment cases. These errors are mainly caused by the interelement incompatibilities, and can

Table 8 Normalized center deflection wc/wref and moments Mc/Mref of SS1 circular plates. Mesh N

3

12

(a) h/R ¼0.02 (h ¼ 0.1) MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9146 0.9573 0.9573 1.0145 1.0790 1.0301

0.9801 0.9900 0.9897 1.0058 1.0075 1.0071

MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9173 1.0453 1.0454 1.0413 0.9646 1.0158

0.9891 1.0085 1.0088 1.0211 0.9995 1.0051

(b) h/R ¼ 0.2 (h ¼ 1) MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9166 0.9535 0.9544 0.9972 1.0766 1.0246

0.9801 0.9880 0.9881 0.9963 1.0096 1.0055

MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9270 1.0473 1.0613 1.0495 0.9787 1.0301

0.9872 1.0143 1.0149 1.0246 1.0062 1.0051

48 wc/wref 0.9951 0.9976 0.9974 1.0014 1.0018 1.0017 Mc/Mref 0.9969 1.0027 1.0030 1.0056 1.0014 1.0015 wc/wref 0.9951 0.9969 0.9969 0.9983 1.0021 1.0014 Mc/Mref 0.9969 1.0046 1.0040 1.0072 1.0010 1.0011

192

Analytical

— — 0.9993 1.0002 1.0043 1.0004

1.0000 (the reference value is 39831.5)

— — 1.0008 1.0014 1.0004 1.0004

1.0000 (the reference value is 5.15625)

— — 0.9992 0.9995 1.0005 1.0004

1.0000 (the reference value is 41.5994)

— — 1.0010 1.0019 1.0003 1.0005

1.0000 (the reference value is 5.15625)

192

Analytical

— — 1.0018 0.9959 0.9968 0.9970

1.0000 (the reference value is 9783.48)

— — 1.0045 1.0012 0.9985 0.9986

1.0000 (the reference value is 2.03125)

— — 1.0010 0.9941 0.9977 0.9971

1.0000 (the reference value is 11.5513)

— — 1.0049 1.0022 0.9983 0.9988

1.0000 (the reference value is 2.03125)

Table 9 Normalized center deflection wc/wref and moments Mc/Mref of clamped circular plates. Mesh N

3

12

(a) h/R ¼0.02 (h ¼ 0.1) MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9269 1.1009 1.1007 0.7604 1.0223 0.8225

0.9914 1.0316 1.0301 0.9455 0.9520 0.9510

MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9255 1.2505 1.2533 0.9607 0.7639 0.8884

1.0092 1.0585 1.0584 1.0152 0.9602 0.9743

(b) h/R ¼ 0.2 (h ¼ 1) MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9311 1.0651 1.0686 0.7359 1.0217 0.8345

0.9897 1.0182 1.0179 0.9205 0.9665 0.9537

MITC4 [2,27] DKMQ [45] ARS-Q12 [14] AC-MQ4 [37] GCP4-12α GCP4-14α

0.9502 1.2603 1.2932 0.9756 0.7972 0.9261

1.0043 1.0732 1.0738 1.0237 0.9762 0.9738

48 wc/wref 0.9983 1.0083 1.0077 0.9859 0.9874 0.9873 Mc/Mref 1.0043 1.0142 1.0166 1.0044 0.9936 0.9940 wc/wref 0.9978 1.0043 1.0041 0.9773 0.9907 0.9883 Mc/Mref 1.0043 1.0191 1.0192 1.0083 0.9925 0.9928

36

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

3.2. The square plate

The square plate is subjected to a uniformly distributed load, as shown in Fig. 4. Due to symmetry, only a quarter of the plate is modeled. The edge length and thickness are denoted by L and h, respectively. Poisson’s ratio μ is set as 0.3. Three cases with different boundary conditions are considered: the clamped case (w¼0, ψn ¼0, ψs ¼0), the soft simply supported (SS1) case (w¼0), and the hard simply supported (SS2) case (w¼ 0, ψs ¼ 0). Tables 2–7 give the dimensionless central deflections and moments of the plate for two ratios (h/L¼ 0.001, 0.1), which are calculated by using the regular mesh (see Fig. 4a). For comparison, results obtained by some other quadrilateral elements [2,10,14,27,37] are also presented. To visually show the convergences, the corresponding plots using log scale are given in Figs. 5–7. Furthermore, the results calculated by using the distorted meshes (see Fig. 4b) are also plotted in Fig. 8. In these figures, the relative errors are defined as follows:     w  wref  M M ref   ; eM ¼   : ew ¼  ð53Þ M ref  wref  Note, because the number of the significant digits for the reference solutions is only four, the low bound of the relative error is about 10  5.

It can be seen that, in most cases, these two elements can provide good results, especially for thin plate cases. But for the SS1 thick plate case, they may converge to a slightly stiffer solution. This problem is mainly caused by the interelement incompatibility. Since the errors are only 1%, the results are still in an acceptable range for engineering applications. An interesting feature is that the results obtained by these two new elements are entirely in agreement with each other in the regular meshes. 3.3. The test for sensitivity to mesh distortions 3.3.1. The clamped thin square plate Fig. 9 shows a clamped thin square plate (h/L¼0.001) subjected to a uniformly distributed load. Due to symmetry, only a quarter of the plate is modeled by a coarse mesh (2  2). Two distortion cases are considered: the symmetric case and the asymmetric case. In the symmetric case, the central mesh node is moved along the main diagonal of the plate to the corner node. In the asymmetric case, the central mesh node is moved to a plate edge along the direction parallel to another plate edge. The normalized results obtained by the new elements and some other elements [2,27,39–41] are plotted in Fig. 10. The results show that the proposed elements are quite insensitive to the mesh distortion, and the element GCP4-14α performs more stably. It

Fig. 16. Distributions of resultants along the radius of clamped circular plate (48 elements).

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

37

should be noted that, in the symmetric case, when the absolute value of the distortion parameter Δ reaches or exceeds 1.25, one element will degenerate into a triangle or a concave quadrangle. Under these situations, the new elements can still perform well, while most other elements fail to provide a solution.

about 1%, those results are still in an acceptable range for engineering applications.

3.3.2. The clamped cantilever plate As shown in Fig. 11, a clamed cantilever plate [42], which is subjected to a tip load P¼ 0.5, is modeled by two elements. The material parameters are E¼ 10.92 and μ¼0.3. The thickness and width are h¼ 0.01 and L¼ 1, respectively. The mesh distortion is measured by the parameter Δ. Fig. 12 shows the percentage deviations of the transverse tip displacement at node C obtained by different elements [14,40,42]. This test proves again that the two new elements are more insensitive to the mesh distortion than others. These two tests sufficiently show the robustness of the new elements in distorted meshes. It can also be found that the element GCP4-14α exhibits a more outstanding performance, since its trial functions are derived from the displacement function with up to the fourth order completeness.

This work is financially supported by the National Natural Science Foundation of China (11272181), the Specialized Research Fund for the Doctoral Program of Higher Education of China (20120002110080), and the Tsinghua University Initiative Scientific Research Program. The authors would also like to thank the support from British Council and the Chinese Scholarship Council through the Sino-UK Higher Education Research Partnership for PhD Studies Scheme, and the support from European Commission through the International Research Staff Exchange Scheme (IRSES).

3.4. The 601 skew plate The Razzaque 601 skew plate [47] (h/L ¼0.001) is subjected to a uniformly distributed load, as shown in Fig. 13. The transverse deflection w and the bending moment My at the central node C are calculated. Results obtained by the new elements and others [2,14,15,27,43–45] are plotted in Fig. 14. The results show that the new elements can produce good results for both deflections and bending moments. 3.5. The circle plate The circular plate subjected to a uniformly distributed load q¼ 1 is shown in Fig. 15. Due to symmetry, only a quarter of the plate is modeled. The geometrical parameters and material parameters are also given in Fig. 15. Two different thickness–radius ratio cases (h/R¼0.02, 0.2) and two different boundary condition cases (the SSl BC case and the clamped BC case) are analyzed. The normalized center deflection and moments are listed in Tables 8 and 9. Fig. 16 shows the distributions of bending moments Mr, Mθ and shear force Tr along the radius of the clamped circular plate calculated by using 48 elements. It can be observed in Fig. 16 that the two new elements can provide resultants which are in good agreement with the exact solutions [43,44], while most conventional displacementbased models fail.

4. Conclusions Two generalized conforming plate elements GCP4-12α and GCP414α, whose trial functions are derived from the analytical solutions of the displacement function [29,30,48], are proposed in this paper. The generalized conforming technique relaxes the requirements of compatibility while improving the element performances. Since the displacement and strain fields of the new elements are assumed according to the fundamental analytical solutions of the displacement function F, which satisfy a priori all related governing equations, these new elements possess the advantages from analytical methods. Some benchmarks are employed to assess the performance of these new elements. In most cases, these two new elements can both provide rather good results, especially for thin plate cases. But for the thick plate under SS1 boundary condition, the result may converge to a slightly stiffer solution. This phenomenon is mainly caused by the interelement incompatibility. Since the error is only

Acknowledgement

References [1] Y.Q. Long, S. Cen, Z.F. Long, Advanced Finite Element Method in Structural Engineering, Springer-Verlag GmbH, Tsinghua University Press, Berlin, Heidelberg, Beijing, 2009. [2] 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. [3] T.J.R. Hughes, R.L. Taylor, W. Kanoknukulchai, A simple and efficient finite element for plate bending, Int. J. Numer. Methods Eng. 11 (10) (1977) 1529–1543. [4] T. Belytschko, C.S. Tsay, W.K. Liu, A stabilization matrix for the bilinear Mindlin plate element, Comput. Meth. Appl. Mech. Eng. 29 (3) (1981) 313–327. [5] W.J. Chen, J.Z. Wang, J. Zhao, Functions for patch test in finite element analysis of the Mindlin plate and the thin cylindrical shell, Sci. China-Phys. Mech. Astron. 52 (5) (2009) 762–767. [6] 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. [7] E. Hinton, H.C. Huang, A family of quadrilateral Mindlin plate element with substitute shear strain fields, Comput. Struct. 23 (3) (1986) 409–431. [8] J.L. Batoz, P. Lardeur, A discrete shear triangular nine dof element for the analysis of thick to very thin plate, Int. J. Numer. Methods Eng. 28 (1989) 533–560. [9] J. Jirousek, A. Venkatesh, Hybrid Trefftz plane elasticity elements with p-method capabilities, Int. J. Numer. Methods Eng. 35 (1992) 1443  1472. [10] O.C. Zienkiewicz, Z.N. Xu, L.F. Zeng, A. Samuelsson, Linked interpolation for Ressiner–Mindlin plate element: Part I—A simple quadrilateral, Int. J. Numer. Methods Eng. 36 (18) (1993) 3043–3056. [11] W.J. Chen, Y.K. Cheung, Refined triangular element based on Mindlin–Reissner plate theory, Int. J. Numer. Methods Eng. 51 (2001) 1259–1281. [12] A. Ibrahimbegović, Quadrilateral finite elements for analysis of thick and thin plates, Comput. Meth. Appl. Mech. Eng. 110 (1993) 195–209. [13] A.K. Soh, Z.F. Long, S. Cen, A new nine DOF triangular element for analysis of thick and thin plates, Comput. Mech. 24 (5) (1999) 408–417. [14] A.K. Soh, S. Cen, Z.F. Long, Y.Q. Long, A new twelve DOF quadrilateral element for analysis of thick and thin plates, Eur. J. Mech. A-Solids 20 (2) (2001) 299–326. [15] W.J. Chen, Y.K. Cheung, Refined quadrilateral element based on Mindlin– Reissner plate theory, Int. J. Numer. Methods Eng. 47 (1–3) (2000) 605–627. [16] G. Castellazzi, P Krysl., A nine-node displacement-based finite element for Reissner–Mindlin plates based on an improved formulation of the NIPE approach, Finite Elem. Anal. Des. 58 (2012) 31–43. [17] P. Hansbo, D. Heintz, M.G. Larson, A finite element method with discontinuous rotations for the Mindlin–Reissner plate model, Comput. Meth. Appl. Mech. Eng. 200 (5–8) (2011) 638–648. [18] P. Hansbo, M.G. Larson, Locking free quadrilateral continuous/discontinuous finite element methods for the Reissner–Mindlin plate, Comput. Meth. Appl. Mech. Eng. 269 (0) (2014) 381–393. [19] D. Ribaric, G. Jelenic, Higher-order linked interpolation in quadrilateral thick plate finite elements, Finite Elem. Anal. Des. 51 (2012) 67–80. [20] N. Nguyen-Thanh, T. Rabczuk, H. Nguyen-Xuan, S. Bordas, An alternative alpha finite element method with discrete shear gap technique for analysis of isotropic Mindlin–Reissner plates, Finite Elem. Anal. Des. 47 (5) (2011) 519–535. [21] G.R. Liu, T. Nguyen-Thoi, Smoothed Finite Element Methods, CRC Press, Taylor and Francis Group, New York, 2010. [22] G.R. Liu, T. Nguyen-Thoi, K.Y. Dai, K.Y. Lam, Theoretical aspects of the smoothed finite element method (SFEM), Int. J. Numer. Methods Eng. 71 (2007) 902–930. [23] H. Nguyen-Xuan, T. Nguyen-Thoi, A stabilized smoothed finite element method for free vibration analysis of Mindlin–Reissner plates, Commun. Numer. Methods Eng. 25 (2009) 882–906.

38

Y. Shang et al. / Finite Elements in Analysis and Design 99 (2015) 24–38

[24] H. Nguyen-Xuan, G.R. Liu, C. Thai-Hoang, T. Nguyen-Thoi, An edge-based smoothed finite element method with stabilized discrete shear gap technique for analysis of Reissner–Mindlin plates, Comput. Meth. Appl. Mech. Eng. 199 (2009) 471–489. [25] T. Nguyen-Thoi, G.R. Liu, H. Nguyen-Xuan, C. Nguyen-Tran, Adaptive analysis using the node-based smoothed finite element method (NS-FEM), Commun. Numer. Methods Eng. 27 (2) (2009) 198–218. [26] T. Nguyen-Thoi, P. Phung-Van, H. Luong-Van, H. Nguyen-Van, H. Nguyen-Xuan, A cell-based smoothed three-node Mindlin plate element (CS-MIN3) for static and free vibration analyses of plates, Comput. Mech. 51 (1) (2013) 65–81. [27] H. Nguyen-Xuan, T. Rabczuk, S. Bordas, J.F. Debongnie, A smoothed finite element method for plate analysis, Comput. Meth. Appl. Mech. Eng. 197 (13) (2008) 1184–1203. [28] D. Ribarić, G. Jelenić, Distortion-immune nine-node displacement-based quadrilateral thick plate finite elements that satisfy constant-bending patch test, Int. J. Numer. Methods Eng. 98 (7) (2014) 492–517. [29] S. Cen, Y. Shang, C.F. Li, H.G. Li, Hybrid displacement function element method: a simple hybrid-Trefftz stress element method for analysis of Mindlin– Reissner plate, Int. J. Numer. Methods Eng. 98 (3) (2014) 203–234. [30] H.C. Hu, Variational Principle of Theory of Elasticity with Applications, Science Press, Gordon and Breach, Science Publisher, Beijing, 1984. [31] Y.Q. Long, K.G. Xin, Generalized conforming element for bending and buckling analysis of plates, Finite Elem. Anal. Des. 5 (1) (1989) 15–30. [32] Y.Q. Long, Y. Xu, Generalized conforming quadrilateral membrane element with vertex rigid rotational freedom, Comput. Struct. 52 (4) (1994) 749–755. [33] Y.Q. Long, Y. Xu, Generalized conforming triangular membrane element with vertex rigid rotational freedoms, Finite Elem. Anal. Des. 17 (4) (1994) 259–271. [34] Y.Q. Long, X.M. Bu, Z.F. Long, Y. Xu, Generalized conforming plate bending elements using point and line compatibility conditions, Comput. Struct. 54 (4) (1995) 717–723. [35] Z.F. Long, Generalized conforming triangular elements for plate bending, Commun. Numer. Methods Eng. 9 (1) (1993) 53–65.

[36] Z.F. Long, Two generalized conforming plate elements based on SemiLoof constraints, Comput. Struct. 47 (2) (1993) 299–304. [37] S. Cen, Y.Q. Long, Z.H. Yao, S.P. Chiew, Application of the quadrilateral area coordinate method: a new element for Mindlin–Reissner plate, Int. J. Numer. Methods Eng. 66 (1) (2006) 1–45. [38] Y.L. Chen, S. Cen, Z.H. Yao, Y.Q. Long, Z.F. Long, Development of triangular flatshell element using a new thin-thick plate bending element based on SemiLoof constrains, Struct. Eng. Mech. 15 (1) (2003) 83–114. [39] Abaqus 6.9, HTML Documentation, Dassault Systèmes Simulia Corp, Providence, RI, USA, 2009. [40] S.L. Weissman, R.L Taylor, Resultant fields for mixed plate bending elements, Comput. Meth. Appl. Mech. Eng. 79 (3) (1990) 321–355. [41] J.L. Batoz, M. Bentahar, Evaluation of a new quadrilateral thin plate bending element, Int. J. Numer. Methods Eng. 18 (11) (1982) 1655  1677. [42] S. de Miranda, F. Ubertini, A simple hybrid stress element for shear deformable plates, Int. J. Numer. Methods Eng. 65 (6) (2006) 808–833. [43] R. Ayad, G. Dhatt, J.L. Batoz, A new hybrid-mixed variational approach for Reissner–Mindlin plates. The MiSP Model, Int. J. Numer. Methods Eng. 42 (1998) 1149–1179. [44] R. Ayad, A. Rigolot, An improved four-node hybrid-mixed element based upon Mindlin’s plate theory, Int. J. Numer. Methods Eng. 55 (2002) 705–731. [45] I. Katili, A new discrete Kirchhoff–Mindlin element based on Mindlin–Reissner plate theory and assumed shear strain fields—Part II: An extended DKQ element for thick-plate bending analysis, Int. J. Numer. Methods Eng. 36 (11) (1993) 1885–1908. [46] H.X. Zhang, J.S. Kuang, Eight-node Reissner–Mindlin plate element based on boundary interpolation using Timoshenko beam function, Int. J. Numer. Methods Eng. 69 (2007) 1345–1373. [47] A. Razzaque, Program for triangular bending elements with derivative smoothing, Int. J. Numer. Methods Eng, 6 (3) (1973) 333  343. [48] Y. Shang, S. Cen, C.F. Li, J.B. Huang, An effective hybrid displacement function element method for solving the edge effect of Mindlin–Reissner plate, Int. J. Numer. Methods Eng. (2014), http://dx.doi.org/10.1002/nme.4843 in press.