Thick Plate Bending Rectangular Element with Shear Deformation
SummaryThis 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 xaxis
θ_{y}  rotation angle about yaxis
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.
FormulationAt a point within the element, coordinate (X,Y), displacement w, rotations θ_{x}, θ_{y} are all interpolated from their nodal values:
y = ∑ N_{i} x_{i}
x = ∑ N_{i} y_{i}
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.
{x_{1}, x_{2}, x_{3}, x_{4}} = {b/2, b/2, b/2, b/2 }
{y_{1}, y_{2}, y_{3}, y_{4}} = {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} = { w_{1}, w_{2}, w_{3}, w_{4} }
Displacements,
u = z θ_{y}
v = z θ_{x}
w
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}
{ε} = [B]{d}, where
{ε} = {ε_{x}, ε_{y}, γ_{xy}, γ_{yz}, γ_{xz}}
{d} = {w_{1}, θ_{x1}, θ_{y1}, w_{2}, θ_{x2}, θ_{y2}, w_{3}, θ_{x3}, θ_{y3}, w_{4}, θ_{x4}, θ_{y4}}
0  0  zN_{1,x}  0  0  zN_{2,x}  0  0  zN_{3,x}  0  0  zN_{4,x}  
0  zN_{1,y}  0  0  zN_{2,y}  0  0  zN_{3,y}  0  0  zN_{4,y}  0  
[B]=[  0  zN_{1,x}  zN_{1,y}  0  zN_{2,x}  zN_{2,y}  0  zN_{3,x}  zN_{3,y}  0  zN_{4,x}  zN_{4,y}  ] (3) 
N_{1,y}  N_{1}  0  N_{2,y}  N_{2}  0  N_{3,y}  N_{3}  0  N_{4,y}  N_{4}  0  
N_{1,x}  0  N_{1}  N_{2,x}  0  N_{2}  N_{3,x}  0  N_{3}  N_{4,x}  0  N_{4} 
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] = [B_{0}] + z[B_{1}]0  0  0  0  0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  0  0  0  0  [0]  
[B_{0}]=[  0  0  0  0  0  0  0  0  0  0  0  0  ] = [  ] (4)  
N_{1,y}  N_{1}  0  N_{2,y}  N_{2}  0  N_{3,y}  N_{3}  0  N_{4,y}  N_{4}  0  [Bg]  
N_{1,x}  0  N_{1}  N_{2,x}  0  N_{2}  N_{3,x}  0  N_{3}  N_{4,x}  0  N_{4} 
0  0  N_{1,x}  0  0  N_{2,x}  0  0  N_{3,x}  0  0  N_{4,x}  
0  N_{1,y}  0  0  N_{2,y}  0  0  N_{3,y}  0  0  N_{4,y}  0  [Bb]  
[B_{1}]=[  0  N_{1,x}  N_{1,y}  0  N_{2,x}  N_{2,y}  0  N_{3,x}  N_{3,y}  0  N_{4,x}  N_{4,y}  ]=[  ]  (5)  
0  0  0  0  0  0  0  0  0  0  0  0  [0]  
0  0  0  0  0  0  0  0  0  0  0  0 
Stiffness Matrix
Stiffness matrix is
[K] = ∫∫∫[B]^{T}[E][B]dxdydz[K] = ∫∫∫([B_{0}]^{T} + z[B_{1}]^{T})[E]([B_{0}] + z[B_{1}])dxdydz
[K] = ∫∫∫([B_{0}]^{T}[E][B_{0}] + z[B_{1}]^{T}[E][B_{0}] + z[B_{0}]^{T}[E][B_{1}] + z^{2}[B_{1}]^{T}[E][B_{1}])dxdydz
Integrate over z (from t/2 to t/2) and noting that ∫z dz from (t/2 to t/2) = 0
[K] = ∫∫([B_{0}]^{T}[E]t[B_{0}] + [B_{1}]^{T}[E](t^{3}/12)[B_{1}])dxdyEt^{3}  1  ν  0  
[Eb]=  [  ν  1  0  ]  (6) 
12(1ν^{2})  0  0  (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]  = ∫∫  [B_{0}]^{T}[D][B_{0}]  dx  dy  
Eb  0  [0]  
= ∫∫  [[0] [Bg]^{T}] [  ][  ] dxdy  
0  Eg  [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
Moments, M_{x} = ∫ σ_{x} z dz M_{y} = ∫ σ_{y} z dz M_{xy} = ∫ τ_{xy} z dz 
{M} = ∫ {σ} z dz
{M} =  [E] ∫ 

z dz 
=  [E] ∫ 

z dz 
=  [E] 

∫ z^{2} dz 
=  [E] 

{Rx} {Ry} 
t^{3} 12 
The above is rearranged as
{M} = [D][Bs]{R}
where
[D] =  Et^{3} 12(1ν^{2}) 

[Bs] = 

{R} = 

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 (Et^{3}/12). However, as thickness increases, the bending stiffness increases much faster. The rate of increase is proportional to t^{2}). 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