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:


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:


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:


Making the same substitution for $u_{i,-1}$, I get:


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?



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.