Formulation of finite element matrices for a 2d elastic isoparametric 4-node quadrilateral element - Eigenvalue analysis of stiffness matrix (energy modes)

Clear["Global`*"]; Off[General::spell, General::spell1];

Shape functions and derivatives

(Subscript[ϕ, 1])  =  (1/4) (1 + (Subscript[ξ, 1])) (1 + (Subscript[ξ, 2]));
(Subscript[ϕ, 2])  =  (1/4) (1 - (Subscript[ξ, 1])) (1 + (Subscript[ξ, 2]));
(Subscript[ϕ, 3])  =  (1/4) (1 - (Subscript[ξ, 1])) (1 - (Subscript[ξ, 2]));
(Subscript[ϕ, 4])  =  (1/4) (1 + (Subscript[ξ, 1])) (1 - (Subscript[ξ, 2]));
dϕ1ξ1  =  
        D[(Subscript[ϕ, 1]), (Subscript[ξ, 1])];
dϕ1ξ2  =  
        D[(Subscript[ϕ, 1]), (Subscript[ξ, 2])];
dϕ2ξ1  =  
        D[(Subscript[ϕ, 2]), (Subscript[ξ, 1])];
dϕ2ξ2  =  
        D[(Subscript[ϕ, 2]), (Subscript[ξ, 2])];
dϕ3ξ1  =  
        D[(Subscript[ϕ, 3]), (Subscript[ξ, 1])];
dϕ3ξ2  =  
        D[(Subscript[ϕ, 3]), (Subscript[ξ, 2])];
dϕ4ξ1  =  
        D[(Subscript[ϕ, 4]), (Subscript[ξ, 1])];
dϕ4ξ2  =  D[(Subscript[ϕ, 4]), (Subscript[ξ, 2])];

Element displacement interpolation matrix

H = (({{(Subscript[ϕ, 1]), 0, (Subscript[ϕ, 2]), 0, (Subscript[ϕ, 3]), 0, (Subscript[ϕ, 4]), 0}, {0, (Subscript[ϕ, 1]), 0, (Subscript[ϕ, 2]), 0, (Subscript[ϕ, 3]), 0, (Subscript[ϕ, 4])}}));
MatrixForm[
      Transpose[H]]

Element nodal displacement vector and dispacement interpolation

Un  =  Table[(Subscript[U, a, i]), {a, 4}, {i, 2}] // Flatten
ue  =  H . Un // Simplify

Nodal coordinates and element coordinates interpolation (mapping)

Xn  =  
      Table[(Subscript[X, a, i]), {a, 4}, {i, 2}] // 
        Flatten
xe  =  H . Xn

Jacobian

J = (({{(∂_((Subscript[ξ, 1])))xe[[1]], (∂_((Subscript[ξ, 1])))xe[[2]]}, {(∂_((Subscript[ξ, 2])))xe[[1]], (∂_((Subscript[ξ, 2])))xe[[2]]}}));
MatrixForm[
      Simplify[Expand[J]]]

Jinv   =  
        Inverse[J] // Simplify;
detJ  =  
      Det[J] // Simplify
A = (1/detJ)(({{J[[2]][[2]], -J[[1]][[2]], 0, 0}, {0, 0, -J[[2]][[1]], J[[1]][[1]]}, {-J[[2]][[1]], J[[1]][[1]], J[[2]][[2]], -J[[1]][[2]]}}));
G = (({{dϕ1ξ1, 0, dϕ2ξ1, 0, dϕ3ξ1, 0, dϕ4ξ1, 0}, {dϕ1ξ2, 0, dϕ2ξ2, 0, dϕ3ξ2, 0, dϕ4ξ2, 0}, {0, dϕ1ξ1, 0, dϕ2ξ1, 0, dϕ3ξ1, 0, dϕ4ξ1}, {0, dϕ1ξ2, 0, dϕ2ξ2, 0, dϕ3ξ2, 0, dϕ4ξ2}}))

Element strain-operator matrix

B  =  A . G // Simplify;
Dimensions[B]

values  =  {(Subscript[X, 1, 1]) -> 1,  (Subscript[X, 1, 2]) -> 1,  
        (Subscript[X, 2, 1]) -> -2,  (Subscript[X, 2, 2]) -> 1,  
        (Subscript[X, 3, 1]) -> -1,  (Subscript[X, 3, 2]) -> -1,  
        (Subscript[X, 4, 1]) -> 1,  
        (Subscript[X, 4, 2]) -> -1}
J  /.  values  //  Simplify // 
      MatrixForm
B  /.  values // Simplify // MatrixForm

Stiffness Matrix

Elastic moduli (Plane strain)

CC = (EE/((1 - (ν^2))))(({{1, ν, 0}, {ν, 1, 0}, {0, 0, (((1 - ν))/2)}}))

values  =  {(Subscript[X, 1, 1]) -> 1,  (Subscript[X, 1, 2]) -> 1,  
        (Subscript[X, 2, 1]) -> -1,  (Subscript[X, 2, 2]) -> 1,  
        (Subscript[X, 3, 1]) -> -1,  (Subscript[X, 3, 2]) -> -1,  
        (Subscript[X, 4, 1]) -> 1,  
        (Subscript[X, 4, 2]) -> -1}
integrand  =  
        detJ Transpose[B] . CC . B   /.  values // 
          Simplify;
K  =  
      Integrate[integrand , {(Subscript[ξ, 1]), -1, 1}, {(Subscript[ξ, 2]), -1, 1}] // 
        Simplify

Eigenvalue analysis of stiffness matrix

k  = 
      (K/((EE/((1 - (ν^2))))))  /.  ν -> 0.3 // 
        Simplify
evalues  =  Eigenvalues[k]
evectors  =  Eigenvectors[k]

(Subscript[U, 1])[v_]  :=  
      Table[evectors[[v]][[j]],  {j, 1, 8, 
          2}]
 u1[v_]  :=  
      Flatten[{(Subscript[U, 1])[v], (Subscript[U, 1])[v][[1]]}]
(Subscript[U, 2])[v_]  :=  
      Table[evectors[[v]][[j]],  {j, 2, 8, 
          2}]
u2[v_]  :=  
      Flatten[{(Subscript[U, 2])[v], (Subscript[U, 2])[v][[1]]}]
Ut[v_]  :=  {u1[v],  u2[v]}
U[v_]  :=  Transpose[Ut[v]]

(Subscript[X, 1]) = Table[Xn[[j]] /. values, {j, 1, 8, 2}];  
    x1  =  Flatten[{(Subscript[X, 1]), (Subscript[X, 1])[[1]]}]
(Subscript[X, 2]) = Table[Xn[[j]] /. values, {j, 2, 8, 2}];  
    x2  =  Flatten[{(Subscript[X, 2]), (Subscript[X, 2])[[1]]}]
X  =  Transpose[{x1, x2}]

b  =  B /.  values  // Simplify;
b . evectors[[8]] // Simplify // Chop

plots  =  
        Table[Show[
            Graphics[{Thickness[0.01], RGBColor[1, 0, 0], Line[X], 
                RGBColor[0, 0, 1], Line[X + U[j]]},  
              AspectRatio -> Automatic, Axes -> True]], {j, 
            8}];