Solving laplace's equation for an inviscid and incompresible fluid Background
I'm working on a 2D inviscid, incompressible fluid sim using vortex methods (that is, treating vortex as discrete particles), and I'm trying to (numerically) solve the no-through boundary conditions, which simply say that the fluid isn't allowed to pass through boundaries of the fluid (which means either the walls of its container or bodies floating in the fluid).  A couple of different sources I have suggest one way of doing it: construct a scalar field such that its gradient stops the fluid flow at the boundaries.  ie:
$$\nabla \Phi \cdot \mathbf{n} = (\mathbf{u_b} - \mathbf{u_\omega}) \cdot \mathbf{n}$$
where $\Phi$ is the scalar field I need to find, $\mathbf{u_\omega}$ is the velocity field of the fluid without the boundary conditions, $\mathbf{u_b}$ is the velocity of the boundary itself, and $\mathbf{n}$ is the normalized surface normal of the boundary.
Because the flow is incompressible, $\Delta \Phi = 0$ (ie: Laplace's equation).  A source I have ("Vortex Methods, Theory and Practice") says you can solve this using a superposition of Green's functions.  This looks like:
$$\frac{\partial \Phi}{\partial \mathbf{n}}(\mathbf{x_b}) = \pm \frac{1}{2} \sigma(\mathbf{x_b}) + \int_S{[\sigma(\mathbf{x_b'}) G(\mathbf{x_b} - \mathbf{x_b'}) -\mu(\mathbf{x_b'}) \mathbf{n} \cdot \nabla G(\mathbf{x_b} - \mathbf{x_b'})] dS(\mathbf{x_b'})}$$
where $\mathbf{x_b}$ is a point on the boundary, $G(\mathbf{x_b} - \mathbf{x_b'})$ is Green's function, and in 2D is $-\frac{\log(|\mathbf{x_b} - \mathbf{x_b'}|)}{2 \pi}$, and $S$ is the boundary I'm trying to prevent the fluid from going through (ie: in 2D it's a polygon's border).  Also:
$$\mu = \Phi_e - \Phi_i$$
That is, it's the difference in the potential ($\Phi$) on opposite sides of the boundary.  The book refers to this as the "double-layer potential".  And:
$$\sigma = \frac{\partial \Phi_e}{\partial \mathbf{n}} - \frac{\partial \Phi_i}{\partial \mathbf{n}} $$
The book refers to this as the single-layer potential or "source".
Then the book assumes that the potential is continuous, so we can say that $\mu = 0$.  A large term in the integral above cancels, and we're finally left with:
Final equation
$$(\mathbf{u_b} - \mathbf{u_\omega}) \cdot \mathbf{n}(\mathbf{x_b}) = -\frac{1}{2} \sigma(\mathbf{x_b}) + \int_S{\sigma(\mathbf{x_b'}) G(\mathbf{x_b} - \mathbf{x_b'}) dS(\mathbf{x_b'})}$$
Questions
First, if anyone has any background in fluids or potential, is the above description correct?  This isn't math I was ever formally taught and I'm just kind of piecing it together as I go.
Also if I understand correctly, in order to actually solve this I need to "discretize" this integral to turn it in to a system of algebraic equations, form a big matrix, and solve it for the $\sigma$ terms.  One simple way of doing this would be to simply select $n$ evenly spaced points along the border of the polygon (is that right?).
However, I don't understand what the $\sigma$ term actually is.  In the definition it looks like it's the Jacobian of $\Phi$ with respect to the normal, which would make it a row vector.  But in the above final equation, it seems to be treating it as a scalar term.
Also, once I have the $\sigma$ terms, what do I actually do with them?  Ultimately I want to calculate $\nabla \Phi (\mathbf{x})$ for some arbitrary point in space $\mathbf{x}$, but I don't see how to do that once you have the $\sigma$ terms.
I found an article on the double-layer potential with math that looks very familiar, so I'm guessing this is a well understood problem, but I'm certainly not understanding it :).  If someone can walk me through the final pieces of math so I can get to a point where I can actually program this, I'd appreciate it.
 A: My understanding is that $\mathbf u_b-\mathbf u_\omega$ is known, and you are solving for $\Phi$. 
In that case, you have the classical  Neumann problem: find the solution of Laplace's equation  with prescribed normal derivative on the boundary. In a formula, 
$$
 \frac{\partial \Phi}{\partial \mathbf n} =g \quad \text{on } S
\tag1$$
where $\frac{\partial \Phi}{\partial \mathbf n}$ is just another way of writing $\nabla \Phi\cdot \mathbf n$, 
and $g=(\mathbf u_b-\mathbf u_\omega)\cdot \mathbf n$ is the desired normal derivative. 
The necessary condition for existence of solution is $\int_S g=0$, since the integral of normal derivative
is equal to  the integral of Laplacian inside, which is zero. 
The relevant potential-theoretic method for the Neumann problem is the single-layer potential. (The double-layer potential is used for the Dirichlet problem; it dropped out from your formulas because it was not needed in the first
place.) It is often explained with analogy from electrostatics: imagine that $\Phi$ is an electrostatic potential; 
then $\nabla \Phi$ is the electric field. Since $\nabla \cdot (\nabla  \Phi)=0$, there must be no charges inside 
the domain. So we place them on the surface $S$ instead. A charge of amount $\sigma$ placed at point $\mathbf p$ of the surface contributes $\sigma (\mathbf p)\,G(\mathbf x-\mathbf p)$ to the electrostatic potential at point $\mathbf x$.  Therefore,
the total potential of the charges distributed along  $S$ with density $\sigma =\sigma(\mathbf p)$ is 
$$
\Phi(x) = \int_S  \sigma(\mathbf p)\, G(\mathbf x-\mathbf p)\, dS(\mathbf x) 
\tag2$$
Now we want to calculate $\nabla \Phi\cdot \mathbf n$  and equate it to $g$. One has to be careful and 
distinguish between two ways in which the charge at $\mathbf p$ contributes to the electric field: 


*

*it creates a  singularity directly at the point $\mathbf p$, where the potential is infinite and its gradient must be 
understood in some generalized sense. 

*it also contributes $\sigma \nabla_{\mathbf x} G(\mathbf x-\mathbf p)$ at points other than $\mathbf p$


This is why the formula for $\nabla \Phi\cdot \mathbf n$  consists of two terms: 
$$
\nabla \Phi(\mathbf x) \cdot \mathbf n  = \frac12 \sigma(\mathbf x) + \int_S  \nabla_{\mathbf x}G (\mathbf x-\mathbf p)\cdot \mathbf n \,d S(\mathbf p) 
\tag3 $$
Some formulas have different $\pm$ signs, which are a matter of convention. But the lack of 
derivative of $G$ in your   formula for the normal derivative of $\Phi$ appears to be an error. 
Since you work with fluid flow, I offer a fluid interpretation: think of placing tiny faucets and sinks along the surface $S$. Counting sinks as negative faucets, we can interpret $\sigma$ as the number of faucets  per unit of area. 
Theoretically, the reformulation of boundary value problem in terms of boundary potentials is helpful because
one can use the well-developed theory of integral operators. But does it help to actually compute the solution to
some precision? I do not know. My (non-expert) impression is that the variational approach, minimizing 
$\int |\nabla u|^2$  among suitable functions with given normal derivative on the boundary, is more practical.   
But you should definitely ask the experts over at Computational Science.    Fluid dynamics is a more active topic there than here.  You may want to either migrate this question there, or post a revised version that focuses on the computational aspect. 
