2
$\begingroup$

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:

enter image description here

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

Regards.

$\endgroup$
3
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.