Mathematica Why is the Neumann boundary condition not satisfied in this FEA simulation?

[member 428835]

Hi PF!

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

 

[member 428835]

Answer: it did not converge because the origin to the curve $\Gamma$ was not at $(0,0)$. Once correcting for this and changing $\hat n$ everything worked out perfectly.