Computing accurate fluxes with FEM

Question

I have solved Poisson equation on a 3d domain with neumann and dirichlet boundary condition. I get the potential, take the gradient for each element and integrate on a surface of an element, I do this for all elements on a surface enclosing an electrode I am interested in knowing the current through. This does not give very accurate results, upto 50% error. And commercial FEM software gets very accurate fluxes.
Does anyone know how they do it or know a better way to compute flux? I am unable to find any references on this

Answer 1

If averaging the element fluxes, as suggested by @Bill Greene, does not produce accurate enough fluxes, there are ways to improve the accuracy of the fluxes by post-processing. A standard and quite straightforward method, which is mentioned as an exercise in the book by Hughes T.J.R., "The Finite Element Method", is as follows:

If your problem is: \begin{align*} & \nabla^2 u = f \quad \text{in} \quad \Omega \\ & \frac{\partial u}{\partial n} = T \quad \text{on} \quad \Gamma_T && \text {Neumann data} \\ & \frac{\partial u}{\partial n} = T^* \quad \text{on} \quad \Gamma_g && \text {sought flux on the Dirichlet boundary} \\ & u = g \quad \text{on} \quad \Gamma_g && \text {Dirichlet data} \end{align*}

Then you find the unknown fluxes by solving the corresponding weak problem: Find $T^*\in L_2(\Gamma_g)$ such that for all $w\in H^{\frac{1}{2}}(\Gamma_g)$ there holds:

\begin{equation*} \int_{\Gamma_g} wT^* d\Gamma = \int_\Omega w_{,i}u_{,i} d\Omega - \int_{\Omega} wf d\Omega - \int_{\Gamma_T} wT d\Gamma, \end{equation*}

where the comma followed by index $i$ denotes differentiation with respect to the coordinates (the summation convention is implied). Note that the Dirichlet data are not used here.

You proceed to discretize this additional weak problem as you did for the main problem.

Answer 2

If $\nabla u$ is of more interest than $u$ itself, it is reasonable to set $\mathbb v := -\sigma \, \nabla u$ and convert your Poisson problem $$ -\nabla \cdot (\sigma \, \nabla u) = f \quad \text{in }\Omega, \\ u = g_E \quad \text{on }\partial\Omega_E, \\ \sigma \, \nabla u \cdot \hat{\mathbb n} = g_N \quad \text{on }\partial\Omega_N $$ to a mixed one: $$ \frac{1}{\sigma} \, \mathbb v + \nabla u = 0 \quad \text{in }\Omega, \\ \nabla \cdot \mathbb v = f \quad \text{in }\Omega, \\ u = g_E \quad \text{on }\partial\Omega_E, \\ \mathbb v \cdot \hat{\mathbb n} = -g_N \quad \text{on }\partial\Omega_N. $$ (Note that for the mixed problem essential BCs are natural and natural BCs are essential.)

Here you build an approximation for $\mathbb v$ directly since it is one of the unknowns. Hence you do not lose accuracy due to numerical differentiation.

This mixed Poisson formulation is pretty well discussed in literature.

Answer 3

One simple way to improve the accuracy of the fluxes is to average the discontinuous values of the element fluxes at each node. This gives you a continuous, nodal flux function over the surface that can be integrated.

If all the elements connected to a node are (roughly) the same size, simply averaging the element fluxes at the node is sufficiently accurate. If the element sizes are significantly different, the contribution of each element flux to the nodal value should be weighted by the size of the element.

There are more sophisticated methods for improving the values of fluxes at the nodes but the fundamental idea is to "smooth" the discontinuous values across inter-element boundaries. In structural analysis, this is commonly referred to as "stress smoothing."

Of course, the fluxes are going to be less accurate than the fundamental nodal unknowns because they are derivatives of these quantities. If you are using the simple 4-node tetrahedron, the flux is constant in the element. Particularly, when flux calculations are important, it is common to use the 10-node tetrahedron element instead.