2
$\begingroup$

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["NDSolve`FEM`"];
(*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, 
    Region`RegionProperty[RegionBoundary[object1], {x, y}, 
       "FastDescription"][[1]][[2]]], 
   DirichletCondition[u[x, y] == V2, 
    Region`RegionProperty[RegionBoundary[object2], {x, y}, 
       "FastDescription"][[1]][[2]]]};
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]
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – user21
    Feb 2, 2023 at 10:45
  • $\begingroup$ 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 $\endgroup$
    – user21
    Feb 2, 2023 at 10:52
  • $\begingroup$ Did you get your application to work? $\endgroup$
    – user21
    Feb 22, 2023 at 10:48

1 Answer 1

2
$\begingroup$

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}];

enter image description here

enter image description here

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Feb 4, 2023 at 0:43
  • $\begingroup$ 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. $\endgroup$
    – user21
    Feb 6, 2023 at 8:50
  • $\begingroup$ BoundaryUnitNormal is now documented in version 13.3 $\endgroup$
    – user21
    Jul 3, 2023 at 9:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.