Formulation of finite element matrices for a 2d elastic isoparametric 4-node quadrilateral element - Motivation for numerical integration

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])}}))

Element nodal displacement vector and dispacement interpolation

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

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]]}}))

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

Elastic moduli (Plane strain)

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

coordinates  =  {(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   /.  coordinates // 
          Simplify;

int  =  integrand[[5]][[1]] // 
        Chop
f  =  
      int  /.  {EE ->  1,  ν -> 
            0.25}
Plot3D[f, {(Subscript[ξ, 1]), -1, 1}, {(Subscript[ξ, 2]), -1, 1}]

inte[r_, s_] :=  
      int  /.  {EE ->  1,  ν -> 0.25,  (Subscript[ξ, 1]) -> 
              r,  (Subscript[ξ, 2]) -> s} // Chop

    inte[(Subscript[ξ, 1]), (Subscript[ξ, 2])]

    Integrate[
        inte[(Subscript[ξ, 1]), (Subscript[ξ, 2])] , {(Subscript[ξ, 1]), -1, 1}, {(Subscript[ξ, 2]), -1, 
          1}] // Simplify

Define the quadrature rule

values  =  {(Subscript[w, 1]) -> 1, 
          (Subscript[w, 2]) -> 
            1, ξ(Subscript[o, 2]) -> -(1/(Sqrt[3])), ξ(Subscript[o, 1]) -> 
            (1/(Sqrt[3]))};
quadrature[
        f_]  :=  (∑_(p = 1)^2)(∑_(q = 1)^2) 
                  (Subscript[w, p]) (Subscript[w, q]) f[ξ(Subscript[o, p]), ξ(Subscript[o, q])]  /.  
        values

quadrature[inte]