Non-linear boundary conditions in electrostatics producing strange results

Hi Everyone,

I have a question relating to solving an electrostatics problem using Poisson's equation, where the problem has a non-linear Neumann boundary condition.
The problem I am attempting to solve is as follows:

$$ \Omega = [0,400]x[0,100]$$

$-\nabla (\epsilon_r\nabla u)=0$ in $\Omega$

$ u = V_g$ on $ \Gamma_{y=0}$

$ \nabla u \cdot n = -13.295 \textrm{sgn}(u)u^2$ on $ \Gamma_{y=100.0}$

$ \nabla u \cdot n = 0.0$ on $ \Gamma_{x=0.0}$

$ \nabla u \cdot n = 0.0$ on $ \Gamma_{x=400.0}$

Hopefully my LaTeX did not fail me in describing the problem...

To explain why I am interested in this problem: I am trying to replicate results for a parallel plate capacitor from a publication (citation below), where one plate (y = 0) has a fixed voltage and the other (y = 100) acquires a charge density whose normal derivative at the interface is described by $-13.295 \textrm{sgn}(u)u^2$. Thus, a non-linear problem.

I believe I have implemented the above problem in FEniCS properly, using a non-linear Newton routine I found in the tutorials. Using the solution for $u$, it is possible to numerically estimate the carrier density using an equation out of [1] $n(x) = 7.3471 * 10^{13} \textrm{sgn}(u)u^2$.

What I expect to see is a square root dependence of carrier density on $V_g$. Unfortunately, instead I am getting a linear dependence, which I find puzzling. I am fairly new to FEnICS (this is the first real problem I am attempting to solve with it), so I am hoping that someone more experienced might be able to suggest where I could have gone wrong. I expect that my issue lies in how I implemented the boundary condition, but I cannot be sure.

Here is my code:

from dolfin import *
import numpy as np
from matplotlib import pyplot as plt

# Inputs/Constants
x_min = 0.0
x_max = 400.0 # device length
x_points = 50

y_min = 0.0
y_max = 100.0 # t_ox = 100 nm
y_points = 50

eps_r = 3.9 # dielectric constant

Vgs = np.arange(-50,51,10.) # Gate voltages
densities = np.zeros(np.size(Vgs)) # Densities to be calculated

# Sub domain for Dirichlet boundary condition
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], x_min)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], x_max)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], y_min)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], y_max)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Create mesh and define function space
mesh = RectangleMesh(Point(x_min,y_min),Point(x_max,y_max),x_points,y_points)

## Initialize mesh function for interior domains
domains = CellFunction("uint", mesh)
domains.set_all(0)

for idx, Vg in enumerate(Vgs): # calculate density for each gate voltage

    # Initialize mesh function for boundary domains
    boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    boundaries.set_all(0)
    left.mark(boundaries, 1)
    top.mark(boundaries, 2)
    right.mark(boundaries, 3)
    bottom.mark(boundaries, 4)

    V = FunctionSpace(mesh, "CG", 1)

    # Define new measures associated with the interior domains and
    # exterior boundaries
    dx = Measure("dx")[domains]
    ds = Measure("ds")[boundaries]

    # Define boundary condition
    bc = DirichletBC(V,Constant(Vg), boundaries,4) # fix bottom gate voltage

    # Define variational problem
    v = TestFunction(V)
    u = Function(V)
    F = inner(grad(u),eps_r*grad(v))*dx \
                  - Constant(0)*v*dx \
                  - 13.295*sign(u)*u*u*v*ds(2)

    # Compute solution
    solve(F == 0, u, bc,solver_parameters={"newton_solver": {"relative_tolerance": 1e-12}})

    # Extract values at graphene layer
    x_points = np.linspace(x_min,x_max,y_points)
    probe_points = np.zeros((y_points,2)) #set of coordinates
    probe_points[:,0] = x_points # x coordinates
    probe_points[:,1] = y_max # y coordinates

    result = np.zeros(y_points)

    for i in xrange(50): # get results at the specified coordinates
        result[i]=u(probe_points[i,:])

    result = np.average(result)
    density = 7.3471E13*np.sign(result)*result**2
    densities[idx] = density

fig, ax = plt.subplots(facecolor='w')
ax.scatter(Vgs,densities)
ax.set_xlabel('Gate voltage')
ax.set_ylabel('Carrier Density')
plt.show()

I kept the grid pretty coarse, so hopefully it runs fast for anyone who tries it. I've done much finer grids and gotten the same results, so that does not seem to be the issue.

The paper, if anyone is curious:

[1] Liu, M. H. (2013). Gate-induced carrier density modulation in bulk graphene: Theories and electrostatic simulation using Matlab pdetool. Journal of Computational Electronics, 12(2), 188–202. http://doi.org/10.1007/s10825-013-0456-9

I would hugely appreciate anyones help! I have been bashing my head against this problem for a while. I expect my problem is something simple and I am just missing it.

Thanks!

Sam