# Solving Laplace for Velocity Potential in Constricted Channel

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?