I am trying to solve the 2D Laplace equation numerically to give the velocity potential of a fluid flowing in a channel with a constriction:
$u_{xx}+u_{yy}=0$
There is a constriction in the channel at $0\leq x\leq 1$, described on the lower boundary ($y=0$) by $y_l(x)=2tx(1-x)$, where $t << 1$. Similarly, at the upper boundary ($y=1$), the constriction is mirrored: $y_u(x)=1-y_l(x)$.
For $x<0$ and $x>1$, the condition at the boundary is given by $\frac{\partial u}{\partial y}=0$.
Along the lower constriction, the boundary is conditioned by $\frac{\partial u}{\partial y}-(1+\frac{\partial u}{\partial x})\frac{dy_l}{dx}=0$ (and equivalently for the upper boundary with $\frac{dy_u}{dx}$). Note, since $t<<1$, I am approximating the boundary to apply to $y=0$ and $y=1$, not on the actual $y_l$ and $y_u$ boundaries.
To start with, I have discretized the Laplace equation, setting the x and y step size to be equal ($\Delta x=\Delta y=h$) to get a stencil for each interior point of the problem:
$u_{i,j}=\frac{1}{4}(u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1})$.
Then, I discretize the $x<0$ and $x>0$ boundary conditions along $y=0$ using central finite differences:
$\frac{\partial u}{\partial y}\approx \frac{u_{i,1}-u_{i,-1}}{2h}, \frac{\partial u}{\partial y}=0 \implies u_{i,1}-u_{i,-1}=0$
Then subbing in from the stencil to eliminate $u_{i,-1}$ as it's outside the domain: $u_{i,0}=\frac{1}{4}\left(2u_{i,1}+u_{i-1,0}+u_{i+1,0}\right)$ (similarly for $y=1$ to get an expression for $u_{i,N}$).
For $0\leq x\leq 1$, I discretize the given boundary condition using central finite differences:
$\frac{u_{i,1}-u_{i,-1}}{2h}-2t(1-2x)-2t(1-2x)\frac{u_{i+1,0}-u_{i-1,0}}{2h}=0$
Making the same substitution for $u_{i,-1}$, I get:
$u_{i,0}=\frac{1}{4}\left(2u_{i,1}+u_{i+1,0}+u_{i-1,0}-2t(1-2x)\left(u_{i+1,0}-u_{i-1,0}\right)-4ht(1-2x)\right)$
and similarly for $u_{i,N}$ for the $y=1$ boundary. At each end on the x axis, I have the condition that $u_x=0$ as $x\to \pm\infty$.
Are these discretizations correct? When implementing them in Python, I get a strange plot:
Darker red means higher potential $u$, white means $u=0$. The two thick red bars indicate $x=0$ and $x=1$ respectively. I would expect my solution plot to somewhat look like the first image in this post but there seems to be a semicircle of 0 potential and of max potential where I would expect just a semicircle of 0 potential (indicating the constriction).
Note, the x and y axis units are grid point index (i,j) rather than the actual x and y values.
What is going wrong?
