For tetrahedral linear elements shape functions Ni can be assumed as equal area coordinates Li given by the formula
where Ve is a volume of tetrahedron. The following integration formula can be useful
∫ | Lα1Lβ2Lμ3Lν4 dxdydz = (α!β!μ!ν!)/(α + β + μ + ν + 3)! 6Ve |
in order to calculate integrals Kb,eac, K*b,ea and Kb,ea. Shape functions for linear elements are Na = La for a = 1, 2, 3, 4. This gives
Kb,ea,c = | ∫ | (∇ Lb)T La∇ Lc dΩe = 1/(4!6Ve)[bb, cb, db] [bc, cc, dc]T = |
1/(4!6Ve)(bbbc + cbcc + dbdc) |
K*b,ea = | ∫ | (∇ Lb)T∇ La dΩe = 1/(36Ve) (bbba + cbca + dbda) |
Kb,ea = | ∫ | La Lb dΩe = 6Ve/5! | when a ≠ b | Ka,ea = | ∫ | La La dΩe = 12Ve/5!. |
Finally, we have
{ | k+ | ( | ∑e Kb,eac | ) | φc + D+ | ∑e | K*b,ea + 1/(Δt) ∑e Kb,ea | } | n+,a = f+,b, |
{ | k- | ( | ∑e Kb,eac | ) | φc + D- ∑e K*b,ea + 1/(Δt) ∑e Kb,ea | } | n-,a = f-,b, |
ε0εφa ∑eKb,ea = |zi|e ∑e Kb,ea (n+,a - n-,a) |
where a, b, c = 1, …, M |
where fi,b = 1/(Δt)(∑e Kb,ea)ni, an-1 with i = {+, -} and only these nodes a, b and c that participate in the particular element e can give a non-zero contribution to the sums of the general type ∑e Ke.
Let us assume that all values of na are known at a time tn. Then we can solve the third equation in equation obtaining the result
φa = (|zi|e)/(ε0ε) {∑eK*b,ea}-1 (∑e Kb,ea) (n+,a - n-,a) |
where {∑eK*b,ea}-1 denotes elements of the inverse matrix to K*. After substitution the solution for φ to equation we get
ni,a {(ki |zi|e)/(ε0ε {∑eK*b,ec}-1 (∑eKb,ec) (n+,c - n-,c) (∑eKb,eac) |
+ Di ∑e K*b,ea + 1/(Δt) ∑e Kb,ea} - fi,b = 0. |
where i={+,-}