# 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.

-

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.

-
Thanks, that's a lot of help. I didn't post in compsci since, while I'm interested in actually implementing this, I'm more interested in the math in the mean time, and this seemed fairly cross-disciplinary. One last question: where does the $\frac{1}{2}\sigma(\mathbf{x})$ come from? I assume it's some limiting case of the integral you gave for $\Phi(\mathbf{x})$? – Jay Lemmon Mar 9 '13 at 23:25
@JayLemmon This is what I tried to explain between (2) and (3). Basically, $\nabla_x G(x-p)$ is not an ordinary function: it is part function, and part a point mass at $p$. The whole thing is called a distribution. When this distribution is integrated, the point mass creates the first term in (3). One can avoid the language of distributions by handling the integral carefully, and this is what potential theory books usually do. You can find this material in Chapter 8 of this book, among many others. – user53153 Mar 9 '13 at 23:44
Assuming you solve for all the $\sigma$ terms, how do you find $\nabla \Phi$? From my course in electromagnetism in college, I'd guess you'd just do a superposition of Coloumb's law for each charge you solved for, which means an inverse square law. Is that right? – Jay Lemmon Mar 10 '13 at 5:12