I am solving the following mixed boundary elasticity problem with the finite element method on the unit circle, using a triangle mesh with 3 nodes on each element:
$ \frac {\partial}{\partial x_j}\sigma_{ij}(u(x_1,x_2))=0 $, $(x_1,x_2) \in \Omega$
$ u_1(x_1,x_2)=u_1^{an}(x_1,x_2)=\frac{1-\nu}{2G(1+\nu)}\sigma_0x_1$ on $\Gamma$
$\sigma_{2j}(u(x_1,x_2))n_j(x_1,x_2)=t_2^{an}(x_1,x_2)$ on $\Gamma$
where:
$ t_2^{an}(x_1,x_2)=\sigma_{21}^{an}(x_1,x_2)n_1+\sigma_{22}^{an}(x1,x2)n_2(x_1,x_2) $ is the 2nd component of the traction vector $ t^{an}=(t_1^{an}, t_2^{an}) $
$ \sigma_i^{an}(x_1,x_2)=\sigma_0 \delta_{ij} $, $i,j=1,2$
$\sigma_0=1.5 *10^{10} N/m^2$ , $ \nu=0.34 $ , $G=3.35*10^{10}N/m^2$
With the following weak form:
$\int_\Omega \sum_{i.j}^2 \sigma_{ij}(u) \epsilon_{ij}(v)= \int_\Gamma v2*t2 $
I have already implemented it correctly for the pure Dirichlet problem. (so the global stiffness matrix is also good for this problem), but I do something wrong when implementing the Neumann boundary condition and I cannot figure out what. (I still get the right solution for u1 where I have the Dirichlet condition).
This is what I do when implementing the Neumann boundary condition:
I go over each triangular element:
-> I extract the coordinates of each node of the element: (x1,y1), (x2,y2), (x3, y3) (the elements are numbered in anticlockwise sense)
-> if (x1,y1) and (x2,y2) are on the boundary:
- I calculate the outward unit normal: $$ n2=\frac{x_1-x_2}{\sqrt{(y_2-y_1)^2+(x_1-x_2)^2}} $$
- I update the force vector F in the place corresponding to $C_e^12$ and $C_e^22$ ($v_2=\sum_{B=1}^3 N^BC_e^B2$) by adding : $$ \sigma_0*n_2*\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}*0.5 $$
-> if (x1,y1) and (x3,y3) are on the boundary:
- $ n2=\frac{x_3-x_1}{\sqrt{(y_1-y_3)^2+(x_3-x_1)^2}} $
- update with $ + (-1)*\sigma_0*n_2*\sqrt{(x_3-x_1)^2+(y_3-y_1)^2}*0.5 $
-> if (x2,y2) and (x3,y3) are on the boundary:
- $ n2=\frac{x_2-x_3}{\sqrt{(y_3-y_2)^2+(x_2-x_3)^2}} $
- update with: $ +\sigma_0*n_2*\sqrt{(x_3-x_2)^2+(y_3-y_2)^2}*0.5 $
The methodology seems fine. The issue could be due to bugs in your code. Check your code carefully and see if there are any bugs.