Welcome to Blogs @ Andrew Qu
Blog Index
All blogs
Search results

Thick Plate Bending Rectangular Element with Shear Deformation


This blog describes the formulation and implementaton of a rectangular plate bending element for thick plates.

Problem Definition

A rectangular plate element is shown in below picture. The element has 4 nodes. There are 3 degrees of freedom at each node,
w - vertical displacement
θx - rotation angle about x-axis
θy - rotation angle about y-axis

Unlike in thin plate where θx, θy are dependent on lateral displacement w, in thick plate, the rotations are independent vairables. We need to formulate a stiffness matrix and detailed implementation.


At a point within the element, coordinate (X,Y), displacement w, rotations θx, θy are all interpolated from their nodal values:

  y = ∑ Ni xi
  x = ∑ Ni yi

where [N] = [N1 N2 N3 N4]                 -----  (1)
  N1 = (1-ξ)(1-η) / 4
  N2 = (1+ξ)(1-η) / 4
  N3 = (1+ξ)(1+η) / 4
  N4 = (1-ξ)(1+η) / 4
are the shape functions. ξ, η range from -1 to 1.

  {x1, x2, x3, x4} = {-b/2, b/2, b/2, -b/2 }
  {y1, y2, y3, y4} = {-h/2, -h/2, h/2, h/2 }
are the coordinates at the corner nodes.

By substituing the nodal coordinates, e found that

  x = (b/2) ξ           ---- (2)
  y = (h/2) η

Similarly, displacement and rotations are
  w = [N]{W}
  θx = [N]{Rx}
  θy = [N]{Ry}

where {W}, {Rx} and {Ry} are the vertical displacement and rotations at corner nodes. For example, {W} = { w1, w2, w3, w4 }

  u = z θy
  v = -z θx


Strains are {ε} = { εx εy γxy γyz γzx}
εx = ∂u/∂x
εy = ∂v/∂y
γxy = ∂u/∂y + ∂v/∂x
γyz = ∂v/∂z + ∂w/∂y
γxz = ∂w/∂x + ∂u/∂z

We ignore εz

εx = ∂u/∂x = z ∂θy/∂x = z ∂[N]/∂x {Ry}
εx = z [N,x]{Ry}
εy = -z [N,y]{Rx}
γxy = z ([N,y]{Ry} - [N,x]{Rx}
γyz = -[N]{Rx} + [N,y]{W}
γxz = [N,x]{W} + [N]{Ry}

Rerrange the terms so that
{ε} = [B]{d}, where
{ε} = {εx, εy, γxy, γyz, γxz}
{d} = {w1, θx1, θy1, w2, θx2, θy2, w3, θx3, θy3, w4, θx4, θy4}

00zN1,x 00zN2,x 00zN3,x 00zN4,x
0-zN1,y0 0-zN2,y0 0-zN3,y0 0-zN4,y0
[B]=[ 0-zN1,xzN1,y 0-zN2,xzN2,y 0-zN3,xzN3,y 0-zN4,xzN4,y] ---(3)
N1,y-N10 N2,y-N20 N3,y-N30 N4,y-N40
N1,x0N1 N2,x0N2 N3,x0N3 N4,x0N4

We will separate matrix [B] into 2 parts, one linear in z and the one containing others. This will be useful later in calculating the stiffness matrix [K].

[B] = [B0] + z[B1]
000 000 000 000
000 000 000 000[0]
[B0]=[ 000 000 000 000] = [     ]  ---(4) 
N1,y-N10 N2,y-N20 N3,y-N30 N4,y-N40[Bg]
N1,x0N1 N2,x0N2 N3,x0N3 N4,x0N4

00N1,x 00N2,x 00N3,x 00N4,x
0-N1,y0 0-N2,y0 0-N3,y0 0-N4,y0[Bb]
[B1]=[ 0-N1,xN1,y 0-N2,xN2,y 0-N3,xN3,y 0-N4,xN4,y]=[ ]  --- (5)
000 000 000 000[0]
000 000 000 000

Stiffness Matrix

Stiffness matrix is

[K] = ∫∫∫[B]T[E][B]dxdydz
[K] = ∫∫∫([B0]T + z[B1]T)[E]([B0] + z[B1])dxdydz
[K] = ∫∫∫([B0]T[E][B0] + z[B1]T[E][B0] + z[B0]T[E][B1] + z2[B1]T[E][B1])dxdydz

Integrate over z (from -t/2 to t/2) and noting that ∫z dz from (-t/2 to t/2) = 0

[K] = ∫∫([B0]T[E]t[B0] + [B1]T[E](t3/12)[B1])dxdy

[Eb]=--------------[ν10]  --- (6)

 1  0
[Eg]=Gt/1.2 [      ] ----- (7)
 0 1

In [Eg], we have introduced an effective shear factor of 1.2, similar to a rectangular beam section.

[Eb] 0
[D] = [      ]
0 [Eg]

[K] = [Kb]+[Kg]

[Kg] = ∫∫ [B0]T[D][B0] dxdy
 Eb0 [0]
 = ∫∫[[0] [Bg]T] [  ][] dxdy
 0Eg [Bg]

 [Kg] = ∫∫[Bg]T[Eg][Bg] dxdy

Replace [Eg] in (7),

 [Kg] = (Gt/1.2)∫∫[Bg]T[Bg] dxdy   --- (8)

 [Kb] = ∫∫[Bb]T[Eb][Bb] dxdy    --- (9)

For the integration in (8), we use 1 Gauss point. For that in (9), we use 4 Gauss points.

Moment Calculation
Stress   {σ} = [E]{ε} where
[E] = E



Mx = -∫ σx z dz
My = -∫ σy z dz
Mxy = -∫ τxy z dz

{M} = -∫ {σ} z dz
{M} = -[E] ∫
z dz
       = -[E] ∫
  z [N,x]{Ry}
-z [N,y]{Rx}
z([N,y]{Ry} - [N,x]{Rx}
z dz
       =   [E]
([N,x]{Rx} - [N,y]{Ry}
∫ z2 dz
       =   [E]
 [0]       -[N,x]
 [N,y]     [0]
([N,x]   -[N,y]


The above is rearranged as

       {M} = [D][Bs]{R}

[D] = Et3

[Bs] =
{R} =
θx1 θx2 θx3 θx4 θy1 θy2 θy3 θy4

Result Comparison

In the test example, we have a 500x500 square plate. The deflection of at the center of the plate is plotted in the diagram below. When the plate is thin, the deflection with shear dformation is smaller than that without shear deformation. This is because, when the plate is thin, shear stiffness (Gt/1.2) stiffens the plate relative to bending stiffness (Et3/12). However, as thickness increases, the bending stiffness increases much faster. The rate of increase is proportional to t2). Therefore, when the plate is thin, the element is too stiff. When the plate is thick, the effect of shear deformation can be upto 4%.


Ads from Googlee
Dr Li Anchor Profi
Engineering anchorage plate design system
©Andrew Qu, 2014. All rights reserved. Code snippets may be used "AS IS" without any kind of warranty. DIY tips may be followed at your own risk.