Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica FEA Neumann BC not satisfied

  1. Jan 9, 2019 #1

    joshmccraney

    User Avatar
    Gold Member

    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"}]
     
    Last edited by a moderator: Jan 9, 2019
  2. jcsd
  3. Jan 16, 2019 at 1:55 PM #2

    joshmccraney

    User Avatar
    Gold Member

    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.
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?