# Implementing Neumann boundary condition for elasticity problem using the finite element method

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$