Four-noded rectangular element shape matrices

Question

I works on a project where I need to compute a modal analysis of an acoustic cavity. The cavity is rigid which translates the problem to the following equation

$$\frac{\partial^2p}{\partial{}x^2}+\frac{\partial^2p}{\partial{}y^2}+k^2p=0$$

and the boundary conditions

$$\partial_xp(y=0)=0,\quad\partial_xp(y=L_y)=0,\quad\partial_yp(x=0)=0,\quad\partial_yp(x=L_x).$$

I'm trying to solve the problem with the finite element method as especially with four-noded rectangular elements. Thanks to the answer to this question I got the form functions which are

\begin{align*} &N_1=\frac{1}{\Delta{}x\Delta{}y}(x-x_2)(y-y_4)\\ &N_2=-\frac{1}{\Delta{}x\Delta{}y}(x-x_1)(y-y_3),\\ &N_3=\frac{1}{\Delta{}x\Delta{}y}(x-x_4)(y-y_2),\\ &N_4=-\frac{1}{\Delta{}x\Delta{}y}(x-x_3)(y-y_1). \end{align*}

Now I'm trying to compute the stiffness and mass matrices of the problem. To do that I used the following Mathematica script:

N1[x_, y_] := (x - x2) (y - y4)/(x1 - x2)/(y1 - y4)
N2[x_, y_] := (x - x1) (y - y3)/(x2 - x1)/(y2 - y3)
N3[x_, y_] := (x - x4) (y - y2)/(x3 - x4)/(y3 - y4)
N4[x_, y_] := (x - x3) (y - y1)/(x4 - x3)/(y4 - y1)
DxN1[x_, y_] := D[N1[x, y], x]
DxN2[x_, y_] := D[N2[x, y], x]
DxN3[x_, y_] := D[N3[x, y], x]
DxN4[x_, y_] := D[N4[x, y], x]
DyN1[x_, y_] := D[N1[x, y], y]
DyN2[x_, y_] := D[N2[x, y], y]
DyN3[x_, y_] := D[N3[x, y], y]
DyN4[x_, y_] := D[N4[x, y], y]
Kx11[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN1[x, y]*DxN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx12[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN1[x, y]*DxN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx13[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN1[x, y]*DxN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx14[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN1[x, y]*DxN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx21[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN2[x, y]*DxN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx22[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN2[x, y]*DxN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx23[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN2[x, y]*DxN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx24[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN2[x, y]*DxN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx31[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN3[x, y]*DxN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx32[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN3[x, y]*DxN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx33[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN3[x, y]*DxN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx34[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN3[x, y]*DxN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx41[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN4[x, y]*DxN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx42[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN4[x, y]*DxN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx43[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN4[x, y]*DxN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Kx44[x_, y_] := 
 Simplify[Integrate[
   Integrate[DxN4[x, y]*DxN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky11[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN1[x, y]*DyN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky12[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN1[x, y]*DyN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky13[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN1[x, y]*DyN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky14[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN1[x, y]*DyN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky21[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN2[x, y]*DyN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky22[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN2[x, y]*DyN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky23[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN2[x, y]*DyN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky24[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN2[x, y]*DyN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky31[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN3[x, y]*DyN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky32[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN3[x, y]*DyN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky33[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN3[x, y]*DyN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky34[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN3[x, y]*DyN4[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky41[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN4[x, y]*DyN1[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky42[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN4[x, y]*DyN2[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky43[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN4[x, y]*DyN3[x, y], {x, x1, x2}], {y, y1, y2}]]
Ky44[x_, y_] := 
 Simplify[Integrate[
   Integrate[DyN4[x, y]*DyN4[x, y], {x, x1, x2}], {y, y1, y2}]]
{{Kx11[x, y], Kx12[x, y], Kx13[x, y], Kx14[x, y]}, {Kx21[x, y], 
   Kx22[x, y], Kx23[x, y], 
   Kx24[x, y]}, {Kx31[x, y], Kx32[x, y], Kx33[x, y], 
    Kx34[x, y]} {Kx41[x, y], Kx42[x, y], Kx43[x, y], 
    Kx44[x, y]}} // MatrixForm
{{Ky11[x, y], Ky12[x, y], Ky13[x, y], Ky14[x, y]}, {Ky21[x, y], 
   Ky22[x, y], Ky23[x, y], 
   Ky24[x, y]}, {Ky31[x, y], Ky32[x, y], Ky33[x, y], 
    Ky34[x, y]} {Ky41[x, y], Ky42[x, y], Ky43[x, y], 
    Ky44[x, y]}} // MatrixForm

I don't like the output because it is not symmetric as it used to be:

Does anyone have the expression of the stiffness and mass matrices for my problem?

Regards.

Answer 1

This is not the best way to make these calculations.

First map the cell $K$

4---------3
|         |
|    K    |
|         |
1---------2

to unit cell [0,1]x[0,1] by $$ \xi = (x - x_1)/\Delta x, \qquad \eta = (y - y_1)/\Delta y $$ The the shape functions are $$ N_1(\xi,\eta) = (1-\xi)(1-\eta), \qquad N_2(\xi,\eta) = \xi (1-\eta), \qquad \ldots, $$ Then $$ \int_K N_i N_j dx dy = \Delta x \Delta y \int_0^1 \int_0^1 N_i N_j d\xi d\eta $$ $$ \int_K \nabla N_i \cdot \nabla N_j dx dy = \Delta x \Delta y \int_0^1 \int_0^1 \left[ \frac{1}{(\Delta x)^2} (\partial_\xi N_i) (\partial_\xi N_j) + \frac{1}{(\Delta y)^2} (\partial_\eta N_i) (\partial_\eta N_j) \right] d\xi d\eta $$ and the integral on the right is easier to compute, even by hand.

This is all elementary stuff, please read some book on FEM before attempting the solution.