Weak Formulation of Poisson's Equation for Electrostatics with Surface Charge

Question

I am currently attempting to use FEnICS to solve an electrostatic problem with two materials of different permittivity $\varepsilon_1$ and $\varepsilon_2$ forming an interface:

Consider a domain $\Omega_1 \cup \Omega_2 = \Omega $ with $\Gamma_T$, $\Gamma_B$, $\Gamma_R$ and $\Gamma_L$ as the top, bottom, right and left exterior boundaries of a simple rectangular domain and $\Gamma$ as the interface between the two dielectrics, the boundary problem is:

$$ -\nabla^2 \varphi = \frac{\rho+\sigma}{\varepsilon} \\ \varphi = U_0 \ \textrm{on} \ \Gamma_T \\ \varphi = U_1 \ \textrm{on} \ \Gamma_B \\ (\nabla \varphi)\cdot\hat{n} = 0 \ \textrm{on} \ \Gamma_R, \ \Gamma_L $$ and finally, the condition due to the surface charge accumulation at the interface: $$ \hat{n}\cdot(D_1-D_2) = \sigma $$ With $D = \varepsilon E$.

Where $\sigma$ is the surface charge. Particularly, I am struggling with the weak formulation of the final boundary condition, I am able to derive, with $w$ being a test function and the general form of $-\nabla^2\varphi = f$ :$$ \int_{\Omega} \nabla\varphi \cdot\nabla w \ d\Omega - \int_{\partial\Omega}w(\nabla\varphi)\cdot\hat{n} \ ds = \int_{\Omega}fw \ d\Omega $$

As I understand, the second integral on the left shall vanish on Dirichlet boundaries and, of course, where a zero Neumann boundary exists.

However, how can I consider the surface charge condition in the weak form on $\Gamma$?

Answer 1

First, I think that we should write the problem as

$$-\nabla(\varepsilon \nabla \phi) = \rho + \sigma\, .$$

Then, the weak form would be

$$\int\limits_\Omega \varepsilon \nabla \phi \nabla w d\Omega - \int\limits_{\partial\Omega} \varepsilon w\nabla\phi\cdot\hat{n}d\Gamma = \int\limits_\Omega f w d\Omega\, .$$

As you mentioned, the second term vanishes for the boundary conditions that you have and you end up with

$$\int\limits_\Omega \varepsilon \nabla \phi \nabla w d\Omega - \int\limits_{\Gamma} \varepsilon w\nabla\phi\cdot\hat{n}d\Gamma = \int\limits_\Omega f w d\Omega\, .$$

But, your interface condition can be rewritten as

$$(-\varepsilon_1 \nabla\phi + \varepsilon_2 \nabla \phi)\cdot\hat{n} = (\varepsilon_1 - \varepsilon_2)\nabla \phi\cdot\hat{n} = \sigma\, ,$$

or

$$\nabla \phi\cdot\hat{n} = \frac{\sigma}{\varepsilon_1 - \varepsilon_2}\, .$$

Plugging this into your weak form

$$\int\limits_\Omega \varepsilon \nabla \phi \nabla w d\Omega = \int\limits_\Omega f w d\Omega + \int\limits_{\Gamma}\frac{w \varepsilon \sigma}{\varepsilon_1 - \varepsilon_2} d\Gamma \, .$$

And I am not sure what value $\varepsilon$ should have for the interface.

Update 20210609

If I understood @ConvexHull properly, you have the following

$$-\nabla^2 \phi = \frac{\rho + \sigma}{\varepsilon}\, .$$

Then, the weak form would be

$$\int\limits_\Omega \nabla \phi \nabla w d\Omega - \int\limits_{\partial\Omega} w\nabla\phi\cdot\hat{n}d\Gamma = \int\limits_\Omega f w d\Omega\, .$$

As you mentioned, the second term vanishes for the boundary conditions that you have and you end up with

$$\int\limits_\Omega \nabla \phi \nabla w d\Omega - \int\limits_{\Gamma} w\nabla\phi\cdot\hat{n}d\Gamma = \int\limits_\Omega f w d\Omega\, .$$

But, your interface condition can be rewritten as

$$(-\varepsilon_1 \nabla\phi + \varepsilon_2 \nabla \phi)\cdot\hat{n} = (\varepsilon_1 - \varepsilon_2)\nabla \phi\cdot\hat{n} = \sigma\, ,$$

or

$$\nabla \phi\cdot\hat{n} = \frac{\sigma}{\varepsilon_1 - \varepsilon_2}\, .$$

Plugging this into your weak form

$$\int\limits_\Omega \nabla \phi \nabla w d\Omega = \int\limits_\Omega f w d\Omega + \int\limits_{\Gamma}\frac{w \sigma}{\varepsilon_1 - \varepsilon_2} d\Gamma \, .$$