# How to specify derivative boundary conditions for the gradient to be normal to the boundary?

I'm currently modeling an electric field with 2 charges. To do so, I use NDSolveValue to solve a Laplacian with 2 Dirichlet conditions on the voltages of the charges, with which I am successfully able to do. However, I now try to add a rectangular conductor with free charges between them, and that's where I encounter my issue.

To set up a conductor, I subtract an additional object from the region and set up conditions on its borders such that the gradient of the field is always perpendicular to the boundary, (as electric field lines going in and out of the conductor are perpendicular to its surface) i.e. for top and bottom $$u_x(x,y) = 0$$, and for right and left $$u_y(x,y)=0$$. I tried using NeumannValue as in code below, but I encounter an error saying

"NDSolveValue::fembcdepderiv: Derivatives of dependent variables in boundary conditions are not supported with the Finite Element Method in this version of NDSolve."


I've looked through similar problems (link link, link), but the solutions proposed either didn't apply to this problem or didn't work for me.

My question is, how can I set up the Neumann boundary condition for the normal derivative to be equal to the gradient, so that the electric field goes perpendicularly inside the conductor? Alternatively, an existing function a-la Acoustic PDE would also work.

Needs["NDSolveFEM"];
(*Define Boundaries*)
boundaries = {{-10, -10}, {10, 10}};

air = Rectangle @@ boundaries;
object1 = Disk[{-8, 0}, 1];
object2 = Disk[{8, 0}, 1];
conductor = Rectangle[{-2, -4}, {2, 4}];
reg12 = RegionUnion[object1, object2, conductor];
reg = RegionDifference[air, reg12];

mesh = ToElementMesh[reg, Transpose[boundaries],
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]];
mesh["Wireframe"]

eq = Laplacian[u[x, y], {x, y}]; V1 = 1; V2 = -1;
bc = {DirichletCondition[u[x, y] == V1,
RegionRegionProperty[RegionBoundary[object1], {x, y},
"FastDescription"][][]],
DirichletCondition[u[x, y] == V2,
RegionRegionProperty[RegionBoundary[object2], {x, y},
"FastDescription"][][]]};
U = NDSolveValue[
{eq ==
NeumannValue[
Derivative[1, 0][u][x,
y], (-4 <= y <= 4 && (x == 2 || 2 + x == 0))]
+ NeumannValue[
Derivative[0, 1][u][x,
y], (-2 <= x <= 2 && (y == 4 || 4 + y == 0))],
bc},
u, {x, y} \[Element] mesh];

ef = -Grad[U[x, y], {x, y}];


Visualization code if you manage to make it work:

{
DensityPlot[U[x, y], {x, y} \[Element] reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotPoints -> 50,
PlotRange -> Transpose[boundaries], ImageSize -> Medium],
Plot3D[U[x, y], {x, y} \[Element] reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotPoints -> 50, PlotRange -> Transpose[boundaries],
ImageSize -> Medium]
}

StreamDensityPlot[Evaluate[ef], {x, y} \[Element] reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> {x, y}, StreamColorFunction -> None,
StreamStyle -> LightGray, PlotRange -> Transpose[boundaries],
ImageSize -> Medium]

• Look at the NeumannValue reference page in the details section for what NeumannValue does. A relation between NeumannValue and Derivative can be found here. I hope that clarifies the meaning and purpose of NeumannValue. Feb 2 at 10:45
• Otherwise you could use BoundaryUnitNormal. You can find an example here. A BoundaryUnitNormal can be part of a NeumannValue - typically the BUNs are generated but you can put them in manually Feb 2 at 10:52
• Did you get your application to work? Feb 22 at 10:48

Is this what you are looking for:

U = NDSolveValue[{eq ==
NeumannValue[
Indexed[BoundaryUnitNormal[x, y],
1], (-4 <= y <= 4 && (x == 2 || 2 + x == 0))] +
NeumannValue[
Indexed[BoundaryUnitNormal[x, y],
2], (-2 <= x <= 2 && (y == 4 || 4 + y == 0))], bc},
u, {x, y} \[Element] mesh];

ef = -Grad[U[x, y], {x, y}];  • This looks very similar to what I'm looking for! I think I can get from here to what I need. Is there documentation for the BoundaryUnitNormal? I've encountered pages where it is present inside of NeumannValue, but I haven't found one that explained how to use it explicitly. Feb 4 at 0:43
• No, currently there is no documentation for BoundaryUnitNormal, mostly because I am lacking applications, so I'd be very interested in seeing your application when you get it to work! Also, if possible, with some background. Feb 6 at 8:50