# Non-linear flux interface condition - variational formulation

Context: I am working on implementing this paper and I am struggling to come up with a variational formulation for the Butler-Volmer interface conditions.

To simplify my question I consider the following problem: \begin{align} -\Delta u &= f &\text{on } \Omega \\ u &= 0 &\text{on } \partial\Omega \\ \frac{\partial u}{\partial n} &= g(u) &\text{on } \Sigma \end{align}

On the following domain, where $\Omega = \Omega_1 \cup \Omega_2$ and $\Sigma$ denotes the interface between the domains. Now I want to derive a variational formulation for this problem, so I can implement it in the FEM software I use (NGSolve). I researched for keywords including mortar and Nitsche methods, transmission conditions and interface conditions, but I'm unsure about the appropriate method and how to adapt it to my problem.

My approach so far:

I write the equation separated for the two domains: \begin{align} -\Delta u_1 &= f &\text{on } \Omega_1 \\ -\Delta u_2 &= f &\text{on } \Omega_2 \\ u_1 &= 0 &\text{on } \partial\Omega \\ u_1 &= u_2 &\text{on } \Sigma \\ \frac{\partial u_1}{\partial n_1} &= g(u_1) &\text{on } \Sigma \\ \frac{\partial u_2}{\partial n_2} &= -g(u_2) &\text{on } \Sigma \end{align}

Then I derive the variational formulation for the outer domain as: \begin{align} \int_{\Omega_1} \nabla u_1 \cdot \nabla v_1 \, dx - \int_{\Sigma} g(u_1) v_1\, ds &= \int_{\Omega_1} f v_1 \, dx \quad \forall v_1 \in H^1_0(\Omega_1) \end{align}

To get a solution for the inner domain I impose the continuity condition on the interface via a mixed formulation: \begin{align} \int_{\Omega_2} \nabla u_2 \cdot \nabla v_2 \, dx + \int_{\Sigma} \lambda v_2 \, ds &= \int_{\Omega_2} f v_2 \, dx \\ \int_{\Sigma} u_2 \mu \, ds &= \int_{\Sigma} u_1 \mu \, ds \end{align}

Question: Am I on the right track? Are there any papers/books considering this problem and outlining different methods?

In case I will have edge element, which evaluates the normal gradient magnitude at integration points on Σ. Note Σ has to have an orientation such that $\mathbf{g}=g(u)\mathbf{n}$, where n is the orientation of internal edge element. Having orientation/direction of edge elements consistent is essential. With that at hand, you can evaluate vector \begin{equation} \mathbf{f}^\Sigma = \sum_{e=1}^{n^e_\Sigma} \int_\Sigma \boldsymbol\varphi^\textrm{T} \left(\mathbf{n}^e \cdot \mathbf{g}(u)\right) \textrm{d}\Sigma \end{equation} where $\boldsymbol\varphi$ is a vector of base functions. $\mathbf{f}^\Sigma$ is added it to the right-hand side. Summation term is to indicate that you loop over edge elements on the internal surface. You have to have assemble $\mathbf{f}^\Sigma$ on the right positions. If you do not like to implement this from scratch, there are couple codes which can do that.
You need also add terms to your tangent matrix to have quadratic convergence. \begin{equation} \mathbf{K}^\Sigma = \sum_{e=1}^{n^e_\Sigma} \int_\Sigma \boldsymbol\varphi^\textrm{T} \mathbf{n}^e \cdot \frac{\partial \mathbf{g}}{\partial{u}} \frac{\partial u}{\partial \overline{u}} \textrm{d}\Sigma = \sum_{e=1}^{n^e_\Sigma} \int_\Sigma \boldsymbol\varphi^\textrm{T} \left(\mathbf{n}^e \cdot \frac{\partial \mathbf{g}}{\partial{u}} \right) \boldsymbol\varphi \textrm{d}\Sigma \end{equation} where $\overline{u}$ are element edge nodal values and $\mathbf{n}^e$ is element edge unit normal. Field $u$ is evaluated using vector of base functions $\boldsymbol\varphi$ as follows \begin{equation} u = \boldsymbol\varphi \cdot \overline{u} \end{equation}