3
$\begingroup$

I'm solving Laplace equation $\nabla^2 \phi = 0$ with BC's $\phi_x(x=\pm 1) = 0,\, \phi_y(y=-h) = 0$ with a specified BC along the circular arc $x^2 + (-1 + y)^2 = 4$, which I call $\Gamma$ (so the domain $\Omega$ is like a rectangular box with the circular arc $\Gamma$ as the lid).

What I've done is taken the analytic solution $\phi_A = \cos(\pi m (x + 1)/2) \cosh(\pi m (y + h)/2) : m \in \mathbb{N}$ and evaluated its normal derivative on $\Gamma$, and then set this as a Neumann boundary condition for the top lid. Since this does not give a unique solution (no Dirichlet BC) I also give the PDE the BC $\phi = \phi_A(1,\Gamma(x=1))$. However, the numerical solution does not converge to the analytic solution. Code is supplied below, which shows the domain $\Omega$ for clarity.

h = 2;
r = 2;
yp = 1;
m = 1;
Ω = ImplicitRegion[-1 <= x <= 1 && -h <= y <= 2 && ! (x^2 + (y - yp)^2 <= r^2), {x, y}];
Region[Ω]
ϕA = Cos[(π m)/2 (x + 1)] Cosh[(π m)/2 (y + h)];
g = ϕA /. y -> -Sqrt[r^2 - x^2] + yp /. x -> 1;
dnϕA = Grad[ϕA, {x, y}]. ({x, y}/Sqrt[x^2 + y^2]);
dnϕΓA = dnϕA /. y -> -Sqrt[r^2 - x^2] + yp;
op = Laplacian[ϕD[x, y], {x, y}];
ΓNV = NeumannValue[dnϕA, x^2 + (y - yp)^2 == r^2];
Γ = {DirichletCondition[ϕD[x, y] == g, x == 1 && y == -Sqrt[r^2 - 1^2] + yp]};
Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.01];
mesh["Wireframe"]
ϕ = NDSolveValue[{op == ΓNV, Γ}, 
ϕD, {x, y} ∈ mesh, InterpolationOrder -> 3];
Plot3D[{ϕ[x, y], ϕA}, {x, y} ∈ Ω, 
 AxesLabel -> {"x", "y", "ϕ"}]
$\endgroup$
  • 2
    $\begingroup$ I think your Neumann value is wrong: since the circle lid is not centered around (0,0), the unit normal in dn\[Phi]A is not (x,y)/norm(...) but (x,y-1)/norm(...), so dn\[Phi]A = Grad[\[Phi]A, {x, y}].(Normalize[{x, y - 1}]); $\endgroup$ – egwene sedai Jan 14 at 22:22
  • 1
    $\begingroup$ @egwenesedai More precisely dn\[Phi]A = Grad[\[Phi]A, {x, y}].({x, y - yp}/Sqrt[x^2 + (y - yp)^2]); $\endgroup$ – Alex Trounev Jan 15 at 13:50
  • $\begingroup$ Thank you egwene sedai and Alex Trounev. I'm new here and unsure how to upvote your comments, but I agree with you both. (Also, a little unclear if I should be thanking each of you...) $\endgroup$ – Josh McCraney Jan 15 at 21:13
  • 2
    $\begingroup$ @JoshMcCraney: something that AFAIK is very welcome on this site in cases where noone answers but people provide useful information in comments is to answer your own question with the information given by others in comments (which also is a good place to mention their help/input)... $\endgroup$ – Albert Retey Jan 15 at 22:28
3
$\begingroup$

egwene sedai and Alex Trounev posted the solution in the comments. I was using the wrong normal vector to the surface when evaluating the NeumannCondition.

Edit: my normal derivative to the circular arc $\Gamma$ assumed the center was the origin. After recognizing the center was intentionally above the origin, the normal vector should be changed to

({x, y - yp}/Sqrt[x^2 + (y - yp)^2])
$\endgroup$
  • 4
    $\begingroup$ It might be worth posting what exactly helped you out so other people don't have to dig through the comments. It's entirely fine etiquette to repost a comment as an answer with attribution. $\endgroup$ – b3m2a1 Jan 16 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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