Hi PF!(adsbygoogle = window.adsbygoogle || []).push({});

I am solving ##\nabla^2\phi = 0## in Mathematica via NDSolveValue. Rather than waste your time explaining the 2D domain ##\Omega##, check out the code at the bottom of page (copy-paste into Mathematica).

Specifically, I am enforcing a Neumann BC on the curved boundary ##\Gamma## through the NeumannValue function. Since the PDE is the Laplace equation, the documentation implies NeumannValue must specify the normal derivative of ##\phi## to ##\Gamma##, or rather ##\nabla \phi \cdot \hat n|_\Gamma##.

I have an analytic solution in the code below I call ##\phi A##, and I compute its normal derivative along ##\Gamma## which I call ##dn\phi\Gamma A##. This is the value I specify the BC at along the curve ##\Gamma##. However, when I run the notebook, this BC is not satisfied, not even when I refine the mesh (see final plot after running code).

Please check out the code below and let me know what you think. I read the documentation for NeumannValue and I think everything looks good but the final plot is still incorrect. Does anyone see the issue?

h = 2;

r = 2;

yp = 1;

m = 1;

\[CapitalOmega] =

ImplicitRegion[-1 <= x <= 1 && -h <= y <=

2 && ! (x^2 + (y - yp)^2 <= r^2), {x, y}];

Region[\[CapitalOmega]]

\[Phi]A = Cos[(\[Pi] m)/2 (x + 1)] Cosh[(\[Pi] m)/2 (y + h)];

g = \[Phi]A /. y -> -Sqrt[r^2 - x^2] + yp /. x -> 1;

dn\[Phi]A = Grad[\[Phi]A, {x, y}]. ({x, y}/Sqrt[x^2 + y^2]);

dn\[Phi]\[CapitalGamma]A = dn\[Phi]A /. y -> -Sqrt[r^2 - x^2] + yp;

op = Laplacian[\[Phi]D[x, y], {x, y}];

\[CapitalGamma]NV = NeumannValue[dn\[Phi]A, x^2 + (y - yp)^2 == r^2];

\[CapitalGamma] = {DirichletCondition[\[Phi]D[x, y] == g,

x == 1 && y == -Sqrt[r^2 - 1^2] + yp]};

Needs["NDSolve`FEM`"]

mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 0.001];

mesh["Wireframe"]

\[Phi] = NDSolveValue[{op == \[CapitalGamma]NV, \[CapitalGamma]}, \

\[Phi]D, {x, y} \[Element] mesh];

n = {x, y}/Sqrt[x^2 + y^2];

dn\[Phi] = Grad[\[Phi][x, y], {x, y}].n;

dn\[Phi]\[CapitalGamma] = dn\[Phi] /. y -> -Sqrt[r^2 - x^2] + yp;

Plot[{dn\[Phi]\[CapitalGamma], dn\[Phi]\[CapitalGamma]A}, {x, -1, 1},

AxesLabel -> {"x",

"\!\(\*SubscriptBox[\(\[Phi]\), \(n\)]\)\!\(\*SubscriptBox[\(|\), \

\(\[CapitalGamma]\)]\)"}, PlotStyle -> {Solid, Dotted},

PlotLegends -> {"numeric", "specified"}]

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Mathematica FEA Neumann BC not satisfied

Have something to add?

**Physics Forums | Science Articles, Homework Help, Discussion**