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

Thick Plate Bending Rectangular Element with Shear Deformation

Summary

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.

Formulation

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 }

Displacements,
  u = z θy
  v = -z θx
  w

Strain

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

Et31ν0
[Eb]=--------------[ν10]  --- (6)
12(1-ν2)00(1-ν)/2

 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]

Now,
[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)

Similarly,
 [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

1-ν2
1ν0
ν10
00(1-ν)/2

Moments,

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

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

12

The above is rearranged as

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

where
[D] = Et3

12(1-ν2)
1ν0
ν10
00(1-ν)/2
[Bs] =
0000-N1,x-N2,x-N3,x-N4,x
N1,yN2,yN3,yN4,y0000
N1,xN2,xN3,xN4,x-N1,y-N2,y-N3,y-N4,y
{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%.

END


Ads from Googlee
Dr Li Anchor Profi
www.anchorprofi.de
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.