Potential flow

Lets say a flow \(\vec{q} = (u,v)\) is imcompressible and irrotational. This means that \(\vec{q}\) must satisfy the following equations:

\begin{equation}\label{1} u_x + v_y = 0 \quad(incompressible) \end{equation} \begin{equation}\label{2} u_y - v_x = 0 \quad(irrotational) \end{equation}

Once a flow field accepts equations \(\ref{1}\) and \(\ref{2}\), we may think think the flow velocity as the gradient of a velocity potential \(\Phi\) (this is what we called potential flow!):

\begin{equation}\label{3} \vec{q} = \nabla\Phi \end{equation}

Now lets substitute equation \(\ref{3}\) back to equations \(\ref{1}\) and \(\ref{2}\) renders:

\begin{equation} \Phi_{xx} + \Phi_{yy} = 0 \end{equation} \begin{equation} \Phi_{xy} - \Phi_{yx} = 0 \end{equation}

It can be seen that potential flow naturally satisfy the irrotational condition. Also, the incompressible condition enforce the velocity potential \(\Phi\) to satisfy Laplace's equation. Thus, the velocity potential \(\Phi\) is a harmonic equation in the domain.

Green's theorem

Before we go further, you may wonder whats the point of introducing Green's theorem here and how it is connected to BEM. Please let me allow the following statement first:

A function, harmonic and continuously differentiable in a closed-regular region R may be represented as the sum of the potentials of a simple and of a double distribution on the boundary of R.

Now let go back to the core idea behind BEM. It only uses the information on boundaries to solves entire domain. Compared with Finite Volume Method (FVM) or Finite Element Method (FEM), BEM helps us saved a lot computational time in engineering problems. Now to prove the above statement, lets consider two arbitrary functions \(G\) and \(F\) in a domain \(A\) bounded by a boundary \(S\) with its normal vector \(\vec{n}\). We can have the following relations (equation 6 is in 2D, but it can be extended to 3D):

\begin{equation}\label{6} \begin{aligned}[t] \int_S(G\frac{\partial F}{\partial \vec{n}}- F\frac{\partial G}{\partial \vec{n}})dS &= \int_S(G\nabla F- F\nabla G)\cdot\vec{n}dS \\ &= \int\int_A \nabla(G\nabla F- F\nabla G)dA \quad \text{Guass divergence theorem applied}\\ &= \int\int_A (\nabla G\nabla F + G\Delta F - \nabla F\nabla G - F\Delta G)dA \\ &= \int\int_A (G\Delta F- F\Delta G)dA \end{aligned} \end{equation}

Equation \(\ref{6}\) is essentially the Green's 2nd identity. Now consider there is a field point P inside the domain and take the function \(G\) as the following function in 2D:

\begin{equation}\label{7} G_{2D} = \frac{\ln{r}}{2\pi} \quad (G_{3D} = -\frac{1}{4\pi r}) \end{equation}

where r is the distance from a the field point P to a control point Q. It can be seen that the function G satisfies the Laplace's equation everywhere when r != 0. As r approach to zero (an infinitesimal disk centered at point P), we have:

\begin{equation} \begin{aligned}[t] \Delta G|_p &= \lim_{r \to 0} \int\int_{A_p} \Delta GdA_p \\ &= \lim_{r \to 0} \oint_{S_p} \nabla G \cdot \vec{n_p} dS_p \\ &= \lim_{r \to 0} \int_0^{2\pi}\frac{1}{2\pi r}r d\theta\\ &= \lim_{r \to 0} 2\pi r \cdot \frac{1}{2\pi r}\\ &= 1\\ \end{aligned} \end{equation}

Thus, we have just proved that function G is a Green's function which satisfy:

\begin{equation} LG(r,r_p) = \delta(r-r_p) \end{equation}

where \(L\) is linear operator such as \(\Delta\) and \(\delta(r-r_p)\) is a delta function which is 1 when the \(r-r_p=0\) and 0 everywhere else. In fact, the function G is a singularity called source. You can try to deploy a source in a domain to understand how source impact a domain at my singularity playground.

For the domain A, everywhere outside the infinitesimal disk \(A_p\) at point P based on equation \(\ref{6}\) satisfies:

\begin{equation} \begin{aligned}[t] \int\int_{A-A_p}({G\Delta F} - \underbrace{F\Delta G}_{I1}) dV &= \int_S(G\frac{\partial F}{\partial \vec{n}}- F\frac{\partial G}{\partial \vec{n}})dS - \int_{S_p}(\underbrace{G\frac{\partial F}{\partial \vec{n_p}}}_{I_2}- \underbrace{F\frac{\partial G}{\partial \vec{n_p}}}_{I_3})dS_p \\ \end{aligned} \end{equation}

where \(\vec{n_p} = -\vec{n}\).

\begin{equation} \begin{aligned}[t] I_1 = 0 \end{aligned} \end{equation} \begin{equation} \begin{aligned}[t] I_2 &= \lim_{r_p \to 0}\frac{\ln{r_p}}{2\pi}\int_{S_p}\frac{\partial F}{\partial \vec{n_p}} dS_p \\ &= \lim_{r_p \to 0} \frac{\ln{r_p}}{2\pi}\int_{S_p} \nabla F\cdot\vec{n_p}dS_p \\ &= \lim_{r_p \to 0} \frac{\ln{r_p}}{2\pi}\int\int_{A_p} \Delta F dV_p \\ &= \lim_{r_p \to 0} \frac{\ln{r_p}}{2\pi}\Delta F\pi r_p^2 \\ &= 0 \end{aligned} \end{equation} \begin{equation} \begin{aligned}[t] I_3 &= \lim_{r_p \to 0}\int_{S_p}F (LG(r_p)) dS_p \\ &= \lim_{r_p \to 0}\int_{S_p}F \delta(r_p) dS_p \\ &= F|_p \end{aligned} \end{equation}

Rerranging terms in equations 10 to 13 renders:

\begin{equation} F|_p = \frac{1}{2\pi}[\int\int_{A}\Delta F \ln{r}dA - \int_S(\ln{r}\frac{\partial F}{\partial n} - \frac{\partial \ln{r}}{\partial n}F)dS] \end{equation}

Equation 14 is the Fredholm integral of the second kind, which means that an unknown depends on the weighted integral of the same unknown function over certain interval. In a 3D case, we replace \(G_{2D}\) to \(G_{3D}\) as equation 7, which renders the Green's third identity:

\begin{equation} F|_p = -\frac{1}{4\pi}\iiint_{V}\frac{\Delta F}{r}\,dV +\frac{1}{4\pi}\iint_{S}\frac{\partial F}{\partial n}\,\frac{1}{r}\,dS -\frac{1}{4\pi}\iint_{S}F\,\frac{\partial}{\partial n}\!\left(\frac{1}{r}\right)\,dS \end{equation} Essentially, equation 14 or 15 tells us that: if we know \(\frac{\partial F}{\partial n}\) and \(F\) on the boundary of the domain, we can solve \(F\) everywhere in the domain!, and this is the core idea behind BEM. Equations 14 and 15 are very general. A great example of its application is in the potential flow problems, where we replace \(F\) with the velocity potential \(\Phi\): \begin{equation} \Phi|_p = \frac{1}{2\pi}\int_S(\ln{r}\frac{\partial \Phi}{\partial n} - \frac{\partial \ln{r}}{\partial n}\Phi)dS \end{equation} In potential flow problem, equation 16 is the governing equation based on total potential formulation for a non-lifting body. To solve equation 16 numerically, we have to discretize the boundary into pieces. Thus, this numerical method is named Boundary Element Method or Panel Method . A example BEM code on this problem can be found here.

In practical application, the governing equation of BEM may be formulated in other ways such as perturbation potential formulation, veclocity based potential formulation, etc. Moreover, BEM can also be applied in structure analysis, such as analyzing the Saint-venant torsion problem (warping), crack prediction, etc.

If you are interested in this topic, Dr. Spyros A Kinnas covers almost everthing on BEM in his course CE380P Boundary Element Methods at UT Austin.