# 2D formulation for Electrostatic Problems

(Difference between revisions)
 Revision as of 14:57, 30 October 2009 (view source)JMora (Talk | contribs)← Older edit Latest revision as of 10:47, 27 November 2009 (view source)JMora (Talk | contribs) (→Source Vector f(e)) (78 intermediate revisions by one user not shown) Line 8: Line 8: \begin{cases} \begin{cases} \left . V - \bar V = 0 \right |_{\Gamma_{V}}  & in ~ \Gamma_{\varphi} \\ \left . V - \bar V = 0 \right |_{\Gamma_{V}}  & in ~ \Gamma_{\varphi} \\ + \, \\ \left . \hat n \vec{D} - \bar D_n = 0 \right |_{\Gamma_{q}}  & in ~ \Gamma_{q} \\ \left . \hat n \vec{D} - \bar D_n = 0 \right |_{\Gamma_{q}}  & in ~ \Gamma_{q} \\ − \left . \frac{\partial V}{\partial r} \right |_{\Gamma_{\infty}} \approx - \frac{V}{r} & in ~ \Gamma_{\infty} + \, \\ + \left . \displaystyle  \frac{\partial V}{\partial r} \right |_{\Gamma_{\infty}} + \approx \displaystyle - \frac{V}{r^{exp}} & in ~ \Gamma_{\infty} \end{cases} \end{cases} [/itex] [/itex] Line 15: Line 18: can be written as (see the [[General_formulation_for_Electrostatic_Problems | General formulation for Electrostatic Problems]]): can be written as (see the [[General_formulation_for_Electrostatic_Problems | General formulation for Electrostatic Problems]]): + + + ::$+ { + \int_{\Omega} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \mathbf{a} \partial \Omega + + \oint_{\Gamma_{\infty}} \mathbf{N^T} \alpha \mathbf{N} \mathbf{a} \partial \Gamma_{\infty} = + \int_{\Omega} \mathbf{N^T} \rho_S \partial \Omega - + \oint_{\Gamma_q} \mathbf{N^T} \bar D_n \partial \Gamma_q - + \oint_{\Gamma_V} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V + } +$ + + + ::$\mathbf{K} \mathbf{a} \,= \mathbf{f}$ Line 23: Line 40: ::$\mathbf{f}^{(e)}= ::[itex]\mathbf{f}^{(e)}= − \int_{\Omega^{(e)}} \mathbf{N^T} \rho_v \partial \Omega^{(e)} - + \int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} - \oint_{\Gamma_q^{(e)}} \mathbf{N^T} \bar D_n \partial \Gamma_q^{(e)} - \oint_{\Gamma_q^{(e)}} \mathbf{N^T} \bar D_n \partial \Gamma_q^{(e)} - \oint_{\Gamma_V^{(e)}} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V^{(e)} \oint_{\Gamma_V^{(e)}} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V^{(e)} Line 31: Line 48: − with: + with ('''''n''''' is the number of nodes of the element): − ::[itex] V (x,y,z) \cong \hat V (x,y,z) = \sum_{i=0}^n N_i (x,y,z) a_i$ + ::$V (x,y) \cong \hat V (x,y) = \sum_{i=0}^n N_i (x,y) a_i = \mathbf{N}^{(e)} · \mathbf{a}^{(e)}$ − ::$\mathbf{N} = + ::[itex]\mathbf{N^{(e)}} = \begin{bmatrix} \begin{bmatrix} N_1 \\ N_1 \\ Line 46: Line 63: \, \\ \, \\ N_n N_n + \end{bmatrix} + \qquad + \mathbf{a^{(e)}} = + \begin{bmatrix} + a_1 \\ + \, \\ + a_2 \\ + \, \\ + \vdots \\ + \, \\ + a_n + \end{bmatrix} + \qquad + \mathbf{B}= \left [ \mathbf{B_1 B_2 ... B_n} \right ] + \qquad + \mathbf{B_i}= + \begin{bmatrix} + \displaystyle \frac{\partial N_i}{\partial x} \\ + \, \\ + \displaystyle \frac{\partial N_i}{\partial y} + \end{bmatrix} + \qquad + \mathbf{\varepsilon}= + \begin{bmatrix} + \varepsilon_x & 0 \\ + \, \\ + 0 & \varepsilon_y \end{bmatrix} \end{bmatrix}$ [/itex] − ::$\mathbf{B_i}= + + ::[itex]\alpha = \frac{1}{|r-r_0|^{exp}} \qquad with \quad exp=0.5, 1, 2...$ + + + + + + + + == 2D formulation for Triangular Elements == + + + After applying the [[Numerical_Integration#Numerical_Integration_for_Isoparametric_Triangular_Domains | numerical integration for triangular elements]] by using the [[Two-dimensional_Shape_Functions#Natural_Coordinates | natural coordinates]], we obtain: + + + ::$+ \mathbf{N^{(e)}} = + \begin{bmatrix} + N_1 & N_2 & N_3 + \end{bmatrix} + = + \begin{bmatrix} + L_1 & L_2 & L_3 + \end{bmatrix} + = + \begin{bmatrix} + (1-\alpha-\beta) & \alpha & \beta + \end{bmatrix} + \qquad + \mathbf{a^{(e)}} = \begin{bmatrix} \begin{bmatrix} − \frac{\partial N_i}{\partial x} \\ + a_1 \\ \, \\ \, \\ − \frac{\partial N_i}{\partial y} \\ + a_2 \\ \, \\ \, \\ − \frac{\partial N_i}{\partial z} + a_3 \end{bmatrix} \end{bmatrix}$ [/itex] + ::$+ \frac{\partial N_1}{\partial \alpha}=-1 \qquad + \frac{\partial N_2}{\partial \alpha}=1 \qquad + \frac{\partial N_3}{\partial \alpha}=0 \qquad + \frac{\partial N_1}{\partial \beta}=-1 \qquad + \frac{\partial N_2}{\partial \beta}=0 \qquad + \frac{\partial N_3}{\partial \beta}=1 +$ − This is: + ::[[Image:NaturalCoordinates 2.jpg|300px]] − ::$− { − \int_{\Omega} W \left [ \nabla^T \mathbf{\varepsilon} \nabla \hat V + \rho_v \right ] \partial\Omega − + \oint_{\Gamma} \overline{W} \left [\mathbf{n}^T \mathbf{\varepsilon} \nabla \hat V + \bar \mathbf{q} + \alpha V \right ] \partial \Gamma=0 − } −$ + ::$+ x = N_1 x_1 + N_2 x_2 + N_3 x_3 = ( 1 - \alpha - \beta) x_1 + \alpha x_2 + \beta x_3 \, +$ − with $\alpha = \frac{1}{r}$ the infinit condition factor and $\bar \mathbf{q}$ the field produced whe $V \,$ is fixed by $\bar V \,$. + ::$+ y = N_1 y_1 + N_2 y_2 + N_3 y_3 = ( 1 - \alpha - \beta) y_1 + \alpha y_2 + \beta y_3 \, +$ − The weak form of this expression can be obtained using the integration by parts. In addition, if     $\bar W = - W \,$: + ::$\mathbf{J^{(e)}} = + \begin{bmatrix} + \displaystyle \frac{\partial x}{\partial \alpha} & \displaystyle \frac{\partial y}{\partial \alpha} \\ \quad \\ + \displaystyle \frac{\partial x}{\partial \beta} & \displaystyle \frac{\partial y}{\partial \beta} + \end{bmatrix} + = + \begin{bmatrix} + - x_1 + x_2 & - y_1 + y_2 \\ + - x_1 + x_3 & - y_1 + y_3 + \end{bmatrix} + \qquad + \mathbf{|J^{(e)}|} = 2 A^{(e)} +$ − ::$− { − \int_{\Omega} \nabla^T W^T \mathbf{\varepsilon} \nabla \hat V \partial \Omega + − \oint_{\Gamma_{\infty}} W^T \alpha V \partial \Gamma_{\infty} = − \int_{\Omega} W^T \rho_v \partial \Omega - − \oint_{\Gamma_q} W^T \bar D_n \partial \Gamma_q - − \oint_{\Gamma_V} W^T \mathbf{q_n} \partial \Gamma_V − } −$ + ::$\mathbf{B(\alpha,\beta)}=\mathbf{J^{(e)}} \mathbf{B(x,y)} \qquad \mathbf{B(x,y)}= \mathbf{[J^{(e)}]^{-1}} \mathbf{B(\alpha,\beta)}$ − Remembering that: + ::$+ \mathbf{B}= + \begin{bmatrix} + \displaystyle \frac{\partial N_1}{\partial x} & + \displaystyle \frac{\partial N_2}{\partial x} & + \displaystyle \frac{\partial N_3}{\partial x}\\ + \, \\ + \displaystyle \frac{\partial N_1}{\partial y} & + \displaystyle \frac{\partial N_2}{\partial y} & + \displaystyle \frac{\partial N_3}{\partial y} + \end{bmatrix} + = + \frac{1}{|\mathbf{J^{(e)}}|} + \begin{bmatrix} + \displaystyle \frac{\partial y}{\partial \beta} & \displaystyle -\frac{\partial y}{\partial \alpha} \\ + \displaystyle -\frac{\partial x}{\partial \beta} & \displaystyle \frac{\partial x}{\partial \alpha} + \end{bmatrix} + \begin{bmatrix} + \displaystyle \frac{\partial N_1}{\partial \alpha} & + \displaystyle \frac{\partial N_2}{\partial \alpha} & + \displaystyle \frac{\partial N_3}{\partial \alpha}\\ + \, \\ + \displaystyle \frac{\partial N_1}{\partial \beta} & + \displaystyle \frac{\partial N_2}{\partial \beta} & + \displaystyle \frac{\partial N_3}{\partial \beta} + \end{bmatrix} +$ − ::$\hat V (x,y,z) = \sum_{i=0}^n N_i (x,y,z) a_i = \mathbf{N} \mathbf{a}^{(e)}$ + ::$+ \mathbf{B} + = + \frac{1}{2 A^{(e)}} + \begin{bmatrix} + - y_1 + y_3 & - y_2 + y_1 \\ + - x_3 + x_1 & - x_1 + x_2 + \end{bmatrix} + \begin{bmatrix} + -1 & 1 & 0 \\ + -1 & 0 & 1 + \end{bmatrix} + = + \frac{1}{2 A^{(e)}} + \begin{bmatrix} + - y_3 + y_2 & - y_1 + y_3 & - y_2 + y_1 \\ + x_3 - x_2 & - x_3 + x_1 & - x_1 + x_2 + \end{bmatrix} +$ − ::$\nabla \hat V = \nabla \mathbf{N} \mathbf{a}^{(e)} = \mathbf{B} \mathbf{a}^{(e)}$ − is the gradient potential with: + === Stiffness Matrix K(e) === + ::$+ \int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)}= + \int \int_{A^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} d x d y = + \int_0^1 \int_0^{1-\beta} |\mathbf{J^{(e)}}| \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} d \alpha d \beta = +$ − ::$\mathbf{B}= \left [ \mathbf{B_1, B_2 ... B_n} \right ]$ + ::$+ = \qquad \qquad |\mathbf{J^{(e)}}| \sum_{p=1}^{n_p} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} W_p = + |\mathbf{J^{(e)}}| \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \sum_{p=1}^{n_p} W_p = + \frac{|\mathbf{J^{(e)}}|}{2} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} +$ − and: + ::::$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = + \begin{bmatrix} + \displaystyle \frac{\partial N_1}{\partial x} & + \displaystyle \frac{\partial N_1}{\partial y} \\ + \, \\ + \displaystyle \frac{\partial N_2}{\partial x} & + \displaystyle \frac{\partial N_2}{\partial y} \\ + \, \\ + \displaystyle \frac{\partial N_3}{\partial x} & + \displaystyle \frac{\partial N_3}{\partial y} + \end{bmatrix} + \begin{bmatrix} + \displaystyle \varepsilon_x & 0 \\ + \, \\ + 0 & \displaystyle \varepsilon_y + \end{bmatrix} + \begin{bmatrix} + \displaystyle \frac{\partial N_1}{\partial x} & + \displaystyle \frac{\partial N_2}{\partial x} & + \displaystyle \frac{\partial N_3}{\partial x}\\ + \, \\ + \displaystyle \frac{\partial N_1}{\partial y} & + \displaystyle \frac{\partial N_2}{\partial y} & + \displaystyle \frac{\partial N_3}{\partial y} + \end{bmatrix} +$ − ::$\mathbf{B_i}= + ::::[itex]\mathbf{B(x,y)^T} \mathbf{\varepsilon} \mathbf{B(x,y)} = + \mathbf{B(\alpha,\beta)^T} \mathbf{[[J^{(e)}]^{-1}]^T} \mathbf{\varepsilon} \mathbf{[J^{(e)}]^{-1}} \mathbf{B(\alpha,\beta)}$ + + + ::::$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = + \frac{1}{|\mathbf{J^{(e)}}|^2} \begin{bmatrix} \begin{bmatrix} − \frac{\partial N_i}{\partial x} \\ + \displaystyle \frac{\partial N_1}{\partial \alpha} & + \displaystyle \frac{\partial N_1}{\partial \beta} \\ \, \\ \, \\ − \frac{\partial N_i}{\partial y} \\ + \displaystyle \frac{\partial N_2}{\partial \alpha} & + \displaystyle \frac{\partial N_2}{\partial \beta} \\ \, \\ \, \\ − \frac{\partial N_i}{\partial z} + \displaystyle \frac{\partial N_3}{\partial \alpha} & + \displaystyle \frac{\partial N_3}{\partial \beta} + \end{bmatrix} + \begin{bmatrix} + \displaystyle \frac{\partial y}{\partial \beta} & \displaystyle -\frac{\partial x}{\partial \beta} \\ + \displaystyle -\frac{\partial y}{\partial \alpha} & \displaystyle \frac{\partial x}{\partial \alpha} + \end{bmatrix} + \begin{bmatrix} + \displaystyle \varepsilon_x & 0 \\ + \, \\ + 0 & \displaystyle \varepsilon_y + \end{bmatrix} + \begin{bmatrix} + \displaystyle \frac{\partial y}{\partial \beta} & \displaystyle -\frac{\partial y}{\partial \alpha} \\ + \displaystyle -\frac{\partial x}{\partial \beta} & \displaystyle \frac{\partial x}{\partial \alpha} + \end{bmatrix} + \begin{bmatrix} + \displaystyle \frac{\partial N_1}{\partial \alpha} & + \displaystyle \frac{\partial N_2}{\partial \alpha} & + \displaystyle \frac{\partial N_3}{\partial \alpha}\\ + \, \\ + \displaystyle \frac{\partial N_1}{\partial \beta} & + \displaystyle \frac{\partial N_2}{\partial \beta} & + \displaystyle \frac{\partial N_3}{\partial \beta} \end{bmatrix} \end{bmatrix}$ [/itex] − The electric field and electric displacement field can be written as follows: + ::::$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = + \frac{1}{(2 A^{(e)})^2} + \begin{bmatrix} + -1 & -1 \\ + 1 & 0 \\ + 0 & 1 + \end{bmatrix} + \begin{bmatrix} + - y_1 + y_3 & - x_3 + x_1 \\ + - y_2 + y_1 & - x_1 + x_2 + \end{bmatrix} + \begin{bmatrix} + \displaystyle \varepsilon_x & 0 \\ + \, \\ + 0 & \displaystyle \varepsilon_y + \end{bmatrix} + \begin{bmatrix} + - y_1 + y_3 & - y_2 + y_1 \\ + - x_3 + x_1 & - x_1 + x_2 + \end{bmatrix} + \begin{bmatrix} + -1 & 1 & 0 \\ + -1 & 0 & 1 + \end{bmatrix} +$ − ::$\mathbf{q} = - \mathbf{B} \mathbf{a}^{(e)} \qquad \mathbf{q'} = - \mathbf{\varepsilon} \mathbf{B} \mathbf{a}^{(e)}$ + ::::$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = + \frac{1}{(2 A^{(e)})^2} + \begin{bmatrix} + - y_3 + y_2 & x_3 + x_2 \\ + - y_1 + y_3 & - x_3 + x_1 \\ + - y_2 + y_1 & - x_1 + x_2 + \end{bmatrix} + \begin{bmatrix} + \displaystyle \varepsilon_x & 0 \\ + \, \\ + 0 & \displaystyle \varepsilon_y + \end{bmatrix} + \begin{bmatrix} + - y_3 + y_2 & - y_1 + y_3 & - y_2 + y_1 \\ + x_3 + x_2 & - x_3 + x_1 & - x_1 + x_2 + \end{bmatrix} +$ − We will now use the Galerkin Method $W_i(x) \equiv N_i(x) \,$. So, finally, the integral expression ready to create the matricial system of equations is: + ::::$+ \int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)}= + A^{(e)} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = + \frac{1}{4 A^{(e)}} + \begin{bmatrix} + - y_3 + y_2 & x_3 + x_2 \\ + - y_1 + y_3 & - x_3 + x_1 \\ + - y_2 + y_1 & - x_1 + x_2 + \end{bmatrix} + \begin{bmatrix} + \displaystyle \varepsilon_x & 0 \\ + \, \\ + 0 & \displaystyle \varepsilon_y + \end{bmatrix} + \begin{bmatrix} + - y_3 + y_2 & - y_1 + y_3 & - y_2 + y_1 \\ + x_3 + x_2 & - x_3 + x_1 & - x_1 + x_2 + \end{bmatrix} +$ − ::$+ ::[itex]\oint_{\Gamma_{\infty}^{(e)}} \mathbf{N^T} \alpha \mathbf{N} \partial \Gamma_{\infty}^{(e)}$ − { + − \int_{\Omega} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \mathbf{a}  \partial \Omega + + − \oint_{\Gamma_{\infty}} \mathbf{N^T} \alpha \mathbf{N} \mathbf{a} \partial \Gamma_{\infty} = + − \int_{\Omega} \mathbf{N^T} \rho_v \partial \Omega - + − \oint_{\Gamma_q} \mathbf{N^T} \bar D_n \partial \Gamma_q - + − \oint_{\Gamma_V} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V + − } + − [/itex] + + === Source Vector f(e) === + ::$+ \int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} = + \int \int_{A^{(e)}} \mathbf{N^T} \rho_S d x d y = + \int_0^1 \int_0^{1-\beta} |\mathbf{J^{(e)}}| \mathbf{N^T} \rho_S d \alpha d \beta = + |\mathbf{J^{(e)}}| \sum_{p=1}^{n_p} \mathbf{N^T} \rho_S W_p = +$ − ::$\mathbf{K} \mathbf{a} \,= \mathbf{f}$ + :::'''Linear case''' ('''np'''=1 integration point): − Note that '''K''' is a coefficients matrix that depends on the geometrical and physical properties of the problem, '''a''' is the vector with the '''''n''''' unknowns to be obtained and '''f''' is a vector that depends on the source values and boundary conditions. + ::::$N=\left [ \frac{1}{3} \quad \frac{1}{3} \quad \frac{1}{3}\right ] \qquad W_i=\frac{1}{2}\,$ − ::$\mathbf{K}^{(e)}= + ::::[itex] − \int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)} + + \int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} − \oint_{\Gamma_{\infty}^{(e)}} \mathbf{N^T} \alpha \mathbf{N} \partial \Gamma_{\infty}^{(e)} + = 2 A^{(e)} \left [ \frac{1}{6} \quad \frac{1}{6} \quad \frac{1}{6}\right ]^T \rho_S + = A^{(e)} \frac{\rho_S}{3} \left [ 1 \quad 1 \quad 1 \right ]^T$ [/itex] − ::$\mathbf{f}^{(e)}= + − \int_{\Omega^{(e)}} \mathbf{N^T} \rho_v \partial \Omega^{(e)} - + − \oint_{\Gamma_q^{(e)}} \mathbf{N^T} \bar D_n \partial \Gamma_q^{(e)} - + :::'''Quadratic case''' ('''np'''=3 integration points): − \oint_{\Gamma_V^{(e)}} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V^{(e)} + −$ + + ::::$p=1 \qquad N=\left [ \frac{1}{2} \quad \frac{1}{2} \quad 0\right ] \qquad W_1=\frac{1}{6}\,$ + ::::$p=2 \qquad N=\left [ 0 \quad \frac{1}{2} \quad \frac{1}{2}\right ] \qquad W_2=\frac{1}{6}\,$ + ::::$p=3 \qquad N=\left [ \frac{1}{2} \quad 0 \quad \frac{1}{2}\right ] \qquad W_3=\frac{1}{6}\,$ + + + ::::$+ \int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} + = 2 A^{(e)} \left ( \left [ \frac{1}{2} \quad \frac{1}{2} \quad 0 \right ]^T + + \left [ 0 \quad \frac{1}{2} \quad \frac{1}{2} \right ]^T + + \left [ \frac{1}{2} \quad 0 \quad \frac{1}{2} \right ]^T \right ) + \frac{1}{6} \rho_S + = A^{(e)} \frac{\rho_S}{3} \left [ 1 \quad 1 \quad 1 \right ]^T +$ + + + + + + + ::$\oint_{\Gamma_q^{(e)}} \mathbf{N^T} \bar D_n \partial \Gamma_q^{(e)}$ + + ::$\oint_{\Gamma_V^{(e)}} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V^{(e)}$ Line 179: Line 450: [[Category: Electrostatic Application]] [[Category: Electrostatic Application]] − [[Category: Theory]] + [[Category: Electrostatic Theory]]

## Latest revision as of 10:47, 27 November 2009

The 2D Electrostatic Poisson's equation given by the governing PDE and its boundary conditions:

$A(V) = \left[ \frac{\partial}{\partial x}\cdot \left( \varepsilon_{x} \cdot \frac{\partial}{\partial x}\right) + \frac{\partial}{\partial y}\cdot \left(\varepsilon_{y} \cdot \frac{\partial}{\partial y} \right) \right]V(x,y) + \rho_S = 0 ~~ in ~ \Omega$

$B(V) = \begin{cases} \left . V - \bar V = 0 \right |_{\Gamma_{V}} & in ~ \Gamma_{\varphi} \\ \, \\ \left . \hat n \vec{D} - \bar D_n = 0 \right |_{\Gamma_{q}} & in ~ \Gamma_{q} \\ \, \\ \left . \displaystyle \frac{\partial V}{\partial r} \right |_{\Gamma_{\infty}} \approx \displaystyle - \frac{V}{r^{exp}} & in ~ \Gamma_{\infty} \end{cases}$

can be written as (see the General formulation for Electrostatic Problems):

${ \int_{\Omega} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \mathbf{a} \partial \Omega + \oint_{\Gamma_{\infty}} \mathbf{N^T} \alpha \mathbf{N} \mathbf{a} \partial \Gamma_{\infty} = \int_{\Omega} \mathbf{N^T} \rho_S \partial \Omega - \oint_{\Gamma_q} \mathbf{N^T} \bar D_n \partial \Gamma_q - \oint_{\Gamma_V} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V }$

$\mathbf{K} \mathbf{a} \,= \mathbf{f}$

$\mathbf{K}^{(e)}= \int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)} + \oint_{\Gamma_{\infty}^{(e)}} \mathbf{N^T} \alpha \mathbf{N} \partial \Gamma_{\infty}^{(e)}$
$\mathbf{f}^{(e)}= \int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} - \oint_{\Gamma_q^{(e)}} \mathbf{N^T} \bar D_n \partial \Gamma_q^{(e)} - \oint_{\Gamma_V^{(e)}} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V^{(e)}$

with (n is the number of nodes of the element):

$V (x,y) \cong \hat V (x,y) = \sum_{i=0}^n N_i (x,y) a_i = \mathbf{N}^{(e)} · \mathbf{a}^{(e)}$

$\mathbf{N^{(e)}} = \begin{bmatrix} N_1 \\ \, \\ N_2 \\ \, \\ \vdots \\ \, \\ N_n \end{bmatrix} \qquad \mathbf{a^{(e)}} = \begin{bmatrix} a_1 \\ \, \\ a_2 \\ \, \\ \vdots \\ \, \\ a_n \end{bmatrix} \qquad \mathbf{B}= \left [ \mathbf{B_1 B_2 ... B_n} \right ] \qquad \mathbf{B_i}= \begin{bmatrix} \displaystyle \frac{\partial N_i}{\partial x} \\ \, \\ \displaystyle \frac{\partial N_i}{\partial y} \end{bmatrix} \qquad \mathbf{\varepsilon}= \begin{bmatrix} \varepsilon_x & 0 \\ \, \\ 0 & \varepsilon_y \end{bmatrix}$

$\alpha = \frac{1}{|r-r_0|^{exp}} \qquad with \quad exp=0.5, 1, 2...$

## 2D formulation for Triangular Elements

After applying the numerical integration for triangular elements by using the natural coordinates, we obtain:

$\mathbf{N^{(e)}} = \begin{bmatrix} N_1 & N_2 & N_3 \end{bmatrix} = \begin{bmatrix} L_1 & L_2 & L_3 \end{bmatrix} = \begin{bmatrix} (1-\alpha-\beta) & \alpha & \beta \end{bmatrix} \qquad \mathbf{a^{(e)}} = \begin{bmatrix} a_1 \\ \, \\ a_2 \\ \, \\ a_3 \end{bmatrix}$

$\frac{\partial N_1}{\partial \alpha}=-1 \qquad \frac{\partial N_2}{\partial \alpha}=1 \qquad \frac{\partial N_3}{\partial \alpha}=0 \qquad \frac{\partial N_1}{\partial \beta}=-1 \qquad \frac{\partial N_2}{\partial \beta}=0 \qquad \frac{\partial N_3}{\partial \beta}=1$

$x = N_1 x_1 + N_2 x_2 + N_3 x_3 = ( 1 - \alpha - \beta) x_1 + \alpha x_2 + \beta x_3 \,$

$y = N_1 y_1 + N_2 y_2 + N_3 y_3 = ( 1 - \alpha - \beta) y_1 + \alpha y_2 + \beta y_3 \,$

$\mathbf{J^{(e)}} = \begin{bmatrix} \displaystyle \frac{\partial x}{\partial \alpha} & \displaystyle \frac{\partial y}{\partial \alpha} \\ \quad \\ \displaystyle \frac{\partial x}{\partial \beta} & \displaystyle \frac{\partial y}{\partial \beta} \end{bmatrix} = \begin{bmatrix} - x_1 + x_2 & - y_1 + y_2 \\ - x_1 + x_3 & - y_1 + y_3 \end{bmatrix} \qquad \mathbf{|J^{(e)}|} = 2 A^{(e)}$

$\mathbf{B(\alpha,\beta)}=\mathbf{J^{(e)}} \mathbf{B(x,y)} \qquad \mathbf{B(x,y)}= \mathbf{[J^{(e)}]^{-1}} \mathbf{B(\alpha,\beta)}$

$\mathbf{B}= \begin{bmatrix} \displaystyle \frac{\partial N_1}{\partial x} & \displaystyle \frac{\partial N_2}{\partial x} & \displaystyle \frac{\partial N_3}{\partial x}\\ \, \\ \displaystyle \frac{\partial N_1}{\partial y} & \displaystyle \frac{\partial N_2}{\partial y} & \displaystyle \frac{\partial N_3}{\partial y} \end{bmatrix} = \frac{1}{|\mathbf{J^{(e)}}|} \begin{bmatrix} \displaystyle \frac{\partial y}{\partial \beta} & \displaystyle -\frac{\partial y}{\partial \alpha} \\ \displaystyle -\frac{\partial x}{\partial \beta} & \displaystyle \frac{\partial x}{\partial \alpha} \end{bmatrix} \begin{bmatrix} \displaystyle \frac{\partial N_1}{\partial \alpha} & \displaystyle \frac{\partial N_2}{\partial \alpha} & \displaystyle \frac{\partial N_3}{\partial \alpha}\\ \, \\ \displaystyle \frac{\partial N_1}{\partial \beta} & \displaystyle \frac{\partial N_2}{\partial \beta} & \displaystyle \frac{\partial N_3}{\partial \beta} \end{bmatrix}$

$\mathbf{B} = \frac{1}{2 A^{(e)}} \begin{bmatrix} - y_1 + y_3 & - y_2 + y_1 \\ - x_3 + x_1 & - x_1 + x_2 \end{bmatrix} \begin{bmatrix} -1 & 1 & 0 \\ -1 & 0 & 1 \end{bmatrix} = \frac{1}{2 A^{(e)}} \begin{bmatrix} - y_3 + y_2 & - y_1 + y_3 & - y_2 + y_1 \\ x_3 - x_2 & - x_3 + x_1 & - x_1 + x_2 \end{bmatrix}$

### Stiffness Matrix K(e)

$\int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)}= \int \int_{A^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} d x d y = \int_0^1 \int_0^{1-\beta} |\mathbf{J^{(e)}}| \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} d \alpha d \beta =$
$= \qquad \qquad |\mathbf{J^{(e)}}| \sum_{p=1}^{n_p} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} W_p = |\mathbf{J^{(e)}}| \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \sum_{p=1}^{n_p} W_p = \frac{|\mathbf{J^{(e)}}|}{2} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B}$

$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = \begin{bmatrix} \displaystyle \frac{\partial N_1}{\partial x} & \displaystyle \frac{\partial N_1}{\partial y} \\ \, \\ \displaystyle \frac{\partial N_2}{\partial x} & \displaystyle \frac{\partial N_2}{\partial y} \\ \, \\ \displaystyle \frac{\partial N_3}{\partial x} & \displaystyle \frac{\partial N_3}{\partial y} \end{bmatrix} \begin{bmatrix} \displaystyle \varepsilon_x & 0 \\ \, \\ 0 & \displaystyle \varepsilon_y \end{bmatrix} \begin{bmatrix} \displaystyle \frac{\partial N_1}{\partial x} & \displaystyle \frac{\partial N_2}{\partial x} & \displaystyle \frac{\partial N_3}{\partial x}\\ \, \\ \displaystyle \frac{\partial N_1}{\partial y} & \displaystyle \frac{\partial N_2}{\partial y} & \displaystyle \frac{\partial N_3}{\partial y} \end{bmatrix}$

$\mathbf{B(x,y)^T} \mathbf{\varepsilon} \mathbf{B(x,y)} = \mathbf{B(\alpha,\beta)^T} \mathbf{[[J^{(e)}]^{-1}]^T} \mathbf{\varepsilon} \mathbf{[J^{(e)}]^{-1}} \mathbf{B(\alpha,\beta)}$

$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = \frac{1}{|\mathbf{J^{(e)}}|^2} \begin{bmatrix} \displaystyle \frac{\partial N_1}{\partial \alpha} & \displaystyle \frac{\partial N_1}{\partial \beta} \\ \, \\ \displaystyle \frac{\partial N_2}{\partial \alpha} & \displaystyle \frac{\partial N_2}{\partial \beta} \\ \, \\ \displaystyle \frac{\partial N_3}{\partial \alpha} & \displaystyle \frac{\partial N_3}{\partial \beta} \end{bmatrix} \begin{bmatrix} \displaystyle \frac{\partial y}{\partial \beta} & \displaystyle -\frac{\partial x}{\partial \beta} \\ \displaystyle -\frac{\partial y}{\partial \alpha} & \displaystyle \frac{\partial x}{\partial \alpha} \end{bmatrix} \begin{bmatrix} \displaystyle \varepsilon_x & 0 \\ \, \\ 0 & \displaystyle \varepsilon_y \end{bmatrix} \begin{bmatrix} \displaystyle \frac{\partial y}{\partial \beta} & \displaystyle -\frac{\partial y}{\partial \alpha} \\ \displaystyle -\frac{\partial x}{\partial \beta} & \displaystyle \frac{\partial x}{\partial \alpha} \end{bmatrix} \begin{bmatrix} \displaystyle \frac{\partial N_1}{\partial \alpha} & \displaystyle \frac{\partial N_2}{\partial \alpha} & \displaystyle \frac{\partial N_3}{\partial \alpha}\\ \, \\ \displaystyle \frac{\partial N_1}{\partial \beta} & \displaystyle \frac{\partial N_2}{\partial \beta} & \displaystyle \frac{\partial N_3}{\partial \beta} \end{bmatrix}$

$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = \frac{1}{(2 A^{(e)})^2} \begin{bmatrix} -1 & -1 \\ 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} - y_1 + y_3 & - x_3 + x_1 \\ - y_2 + y_1 & - x_1 + x_2 \end{bmatrix} \begin{bmatrix} \displaystyle \varepsilon_x & 0 \\ \, \\ 0 & \displaystyle \varepsilon_y \end{bmatrix} \begin{bmatrix} - y_1 + y_3 & - y_2 + y_1 \\ - x_3 + x_1 & - x_1 + x_2 \end{bmatrix} \begin{bmatrix} -1 & 1 & 0 \\ -1 & 0 & 1 \end{bmatrix}$

$\mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = \frac{1}{(2 A^{(e)})^2} \begin{bmatrix} - y_3 + y_2 & x_3 + x_2 \\ - y_1 + y_3 & - x_3 + x_1 \\ - y_2 + y_1 & - x_1 + x_2 \end{bmatrix} \begin{bmatrix} \displaystyle \varepsilon_x & 0 \\ \, \\ 0 & \displaystyle \varepsilon_y \end{bmatrix} \begin{bmatrix} - y_3 + y_2 & - y_1 + y_3 & - y_2 + y_1 \\ x_3 + x_2 & - x_3 + x_1 & - x_1 + x_2 \end{bmatrix}$

$\int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)}= A^{(e)} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} = \frac{1}{4 A^{(e)}} \begin{bmatrix} - y_3 + y_2 & x_3 + x_2 \\ - y_1 + y_3 & - x_3 + x_1 \\ - y_2 + y_1 & - x_1 + x_2 \end{bmatrix} \begin{bmatrix} \displaystyle \varepsilon_x & 0 \\ \, \\ 0 & \displaystyle \varepsilon_y \end{bmatrix} \begin{bmatrix} - y_3 + y_2 & - y_1 + y_3 & - y_2 + y_1 \\ x_3 + x_2 & - x_3 + x_1 & - x_1 + x_2 \end{bmatrix}$

$\oint_{\Gamma_{\infty}^{(e)}} \mathbf{N^T} \alpha \mathbf{N} \partial \Gamma_{\infty}^{(e)}$

### Source Vector f(e)

$\int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} = \int \int_{A^{(e)}} \mathbf{N^T} \rho_S d x d y = \int_0^1 \int_0^{1-\beta} |\mathbf{J^{(e)}}| \mathbf{N^T} \rho_S d \alpha d \beta = |\mathbf{J^{(e)}}| \sum_{p=1}^{n_p} \mathbf{N^T} \rho_S W_p =$

Linear case (np=1 integration point):

$N=\left [ \frac{1}{3} \quad \frac{1}{3} \quad \frac{1}{3}\right ] \qquad W_i=\frac{1}{2}\,$
$\int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} = 2 A^{(e)} \left [ \frac{1}{6} \quad \frac{1}{6} \quad \frac{1}{6}\right ]^T \rho_S = A^{(e)} \frac{\rho_S}{3} \left [ 1 \quad 1 \quad 1 \right ]^T$

$p=1 \qquad N=\left [ \frac{1}{2} \quad \frac{1}{2} \quad 0\right ] \qquad W_1=\frac{1}{6}\,$
$p=2 \qquad N=\left [ 0 \quad \frac{1}{2} \quad \frac{1}{2}\right ] \qquad W_2=\frac{1}{6}\,$
$p=3 \qquad N=\left [ \frac{1}{2} \quad 0 \quad \frac{1}{2}\right ] \qquad W_3=\frac{1}{6}\,$

$\int_{\Omega^{(e)}} \mathbf{N^T} \rho_S \partial \Omega^{(e)} = 2 A^{(e)} \left ( \left [ \frac{1}{2} \quad \frac{1}{2} \quad 0 \right ]^T + \left [ 0 \quad \frac{1}{2} \quad \frac{1}{2} \right ]^T + \left [ \frac{1}{2} \quad 0 \quad \frac{1}{2} \right ]^T \right ) \frac{1}{6} \rho_S = A^{(e)} \frac{\rho_S}{3} \left [ 1 \quad 1 \quad 1 \right ]^T$

$\oint_{\Gamma_q^{(e)}} \mathbf{N^T} \bar D_n \partial \Gamma_q^{(e)}$
$\oint_{\Gamma_V^{(e)}} \mathbf{n^T} \mathbf{N^T} \mathbf{q_n} \partial \Gamma_V^{(e)}$