# Vortex panel solution method

This section will use vortex panels to solve to the flow around general objects by setting up and solving a system of linear equations.

## General vortex sheets

A curved vortex sheet with a variable strength can describe the flow around any immersed object. This is achieved by having the sheet act as an infinitely thin version of the boundary layer to enforce the no-slip boundary condition. 

---



---

In other words we use the sheets to force the tangential velocity $u_s$ to zero at every point $s$ on the body surface $\cal S$

\begin{equation}
u_s = \vec u \cdot \hat s = 0 \quad s \in \cal S
\end{equation}

Substituting the definition of velocity induced by the sheet, the tangential velocity condition is

\begin{equation}
\left[\vec U+\frac{\partial}{\partial \vec x}\oint_{\cal S} \frac{\gamma(s')}{2\pi}\theta(s,s')\ ds'\right] \cdot\hat s = 0 
\end{equation}

where $\vec U$ is the background velocity that has been added by superposition. 

**If we can use this equation to determine the strength distribution $\gamma(s)$ along the sheet then we will have solved for the potential flow around the body!**

## Discrete vortex panels

For general body surface shapes the velocity is a highly nonlinear function of $\gamma(s)$, rendering analytic solution unlikely. We could attempt some complex analytic expansions, but why would we want to do that?

##### Numerical fundamental: Discretization
##### Replace continuous functions with linear approximations

We already know that the velocity depends **linearly** on $\gamma$ for a vortex panel. This makes it easy to solve for $\gamma$ as a function of $u_s$.

If we break up the continuous sheet into a series of vortex panels, we can add their influence together using superposition. We can then apply the no slip boundary condition, and use it to solve for $\gamma$. 




This is the essence of the *vortex panel method*.

---

## Array of panels

To help make this more concrete, lets use a unit circle as an example. The `make_circle` function has already been implemented in the `VortexPanel` module. Since we're going to use many functions in that module lets import the whole thing (like we do for `numpy`), and since it has a long name, lets give it the nick-name `vp`:


```python
# Annoying stuff to find files in other directories (only needed once)
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)

# Actual import code
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
from vortexpanel import VortexPanel as vp
```

Since you need to build your own shapes for the project, lets `inspect` how `make_circle` works


```python
from inspect import getsource
print(getsource(vp.make_circle))
```

Most of that is a comment block explaining the inputs and giving an example of using the function. Only the last three lines are code. An array of `theta` values are defined and then converted into `x,y` positions. The last line converts those points into an array of vortex panels using a function called `panelize`. 

Lets use the two lines in the example to make a circle of panels and plot it:


```python
circle = vp.make_circle(N=32)  # make
circle.plot(style='o-')        # plot
plt.axis('equal');             # show its a circle
```

Play around with the `style` argument to see the points and panels that make up the shape. 

Note that in addition to holding the array of `panels`, the `PanelArray` object also has some really useful functions like `PanelArray.plot()`.

The `PanelArray` class (and `Panel` within it) is written in an [Object Oriented Programming](https://realpython.com/python3-object-oriented-programming/#what-is-object-oriented-programming-oop) style.

##### Coding fundamental: Objects
##### Object-Oriented Programming helps organize data and functions

Lets see what the help says


```python
help(vp.PanelArray)
```

From the help, we can see that `PanelArray` contains the array of panels as well as an angle off attack. `Panel` also contains functions (the `methods`) which can be applied to it`self` such as plotting the panel and determining the panel's induced velocity. This avoids needing to pass the data about the geometry to the function - it's already all built-in!

---
Now that we have an example, what is the velocity induced by these panels and the uniform flow? 

Using superposition, the total velocity at any point $x,y$ is simply

\begin{equation}
\vec u(x,y) = \vec U+\sum_{j=0}^{N-1} \vec u_j(x,y)
\end{equation}

where we use the index  $j$ to label each of the $N$ panels. This is easily implemented for an array of panels:


```python
print(getsource(vp.PanelArray.velocity))
```

where the `for-in` notation loops through every panel in the array.

---
The function `PanelArray.plot_flow(size=2)` uses this superposition code to visualize the flow:


```python
circle.plot_flow()
```

##### Quiz

Why is the flow going through the body above?

1. We haven't applied the flow tangency condition
1. We haven't applied the no-slip condition
1. We haven't determined the correct `gamma` for each panel

## System of linear equations

One key to the solution method is that the velocity $\vec u_j$ induced by vortex panel $j$ depends __linearly__ on $\gamma_j$. So lets write $\vec u_j$ in a way that makes the linearity explicit:

$$ \vec u_j(x,y)=\gamma_j\ \vec f_j(x,y)$$

where $\vec f_j$ is the part of the velocity function which depends only on the panel geometry.

Substituting this linear relation for the velocity $\vec u_j$ into the total velocity equation and applying the no-slip boundary condition we have:

$$ u_s = \left[\vec U + \sum_{j=0}^{N-1} \gamma_j \ \vec f_j(x,y)\right]\cdot\hat s=0 $$

So the goal is to set $\gamma$ on each panel such that this condition is enforced on the body.

##### Quiz

How many unknowns are there?

1. $1$
1. $N$
1. $N^2$

But we only have one equation, the no-slip condition... right?

##### Numerical fundamental: Consistency
##### Develop enough equations to match the unknowns

We need to have as many equations as we have unknowns to be consistent and to be able to determine a solution.

Luckily the no-slip condition is a continuous equation - it applies to *every* point on the body. Therefore, we can evaluate the boundary equation **on each panel**. Then we will have a consistent linear system.

---

Let's start with one panel, say panel $i$, and set the total tangential velocity at the center of the panel to zero:

$$ \frac 12 \gamma_i + \vec U\cdot\hat s_i + \sum_{j=0, j\ne i}^{N-1} \gamma_j \ \vec f_j(x_i,y_i)\cdot\hat s_i= 0 $$

Note that we've used the simple relation for the velocity a panel induces on itself

$$ \vec u_i(x_i,y_i) \cdot \hat s_i = \frac 12 \gamma_i $$

for the tangential velocity that the panel induces on itself.

Next, let's write the summation as an inner product of two arrays in order to separate out the knowns from the unknowns:

$$
\begin{bmatrix} \vec f_0(x_i,y_i)\cdot\hat s_i & \vec f_1(x_i,y_i)\cdot\hat s_i & \cdots & \frac 12 & \cdots & \vec f_{N-1}(x_i,y_i)\cdot\hat s_i\end{bmatrix} \times \begin{bmatrix} \gamma_0 \\ \gamma_1 \\ \vdots \\ \gamma_i \\ \vdots \\ \gamma_{N-1} \end{bmatrix} + \vec U \cdot \hat s_i = 0
$$

Written like this, we can see two things:

 - The no-slip condition on panel $i$ depends on the strength at every panel, and
 - This is just the $i$th row of a matrix of equations

\begin{equation}
\begin{bmatrix} 
\frac 12  & \vec f_1(x_0,y_0)\cdot\hat s_0 & \cdots & \vec f_{N-1}(x_0,y_0)\cdot\hat s_0\\[0.5em]
\vec f_0(x_1,y_1)\cdot\hat s_1 & \frac 12 & \cdots & \vec f_{N-1}(x_1,y_1)\cdot\hat s_1 \\[0.5em]
\vdots & \vdots & \ddots & \vdots \\[0.5em] 
\vec f_0(x_{N-1},y_{N-1})\cdot\hat s_{N-1} & \vec f_1(x_{N-1},y_{N-1})\cdot\hat s_{N-1} & \cdots & \frac 12
\end{bmatrix} 
\times \begin{bmatrix} \gamma_0 \\[0.9em] \gamma_1 \\[0.9em] \vdots \\[0.9em] \gamma_{N-1} \end{bmatrix} 
= -\begin{bmatrix} \vec U\cdot\hat s_0 \\[0.7em] \vec U\cdot\hat s_1 \\[0.7em] \vdots \\[0.7em] \vec U\cdot\hat s_{N-1} \end{bmatrix}
\end{equation}

which we can summarize as: 

$$ A \gamma = b$$

this defines the complete linear system.

---

## Summary

Lets review the key concepts:

1. The no-slip vortex sheet equation implicitly defines the analytic $\gamma(s)$, but is highly nonlinear.
1. Breaking the sheet up into a set of panels ensures the velocity is a linear function of $\gamma_i$. 
1. Applying the no-slip condition to each panel results in a linear system of equations for $\gamma_i$. 

Therefore, the solution $\gamma_i$ of $A\gamma=b$ is a numerical approximation to the analytic solution $\gamma(s)$. 

*This same basic procedure is applied whenever numerically approximating partial differential equations!*

---

## Implementation

We can construct `A` and `b` in only a few lines of code:
```python
    # construct the matrix `A`
    for j, p_j in enumerate(panels):      # loop through columns
        fx,fy = p_j.constant(xc,yc)           # f_j at all panel centers
        A[:,j] = fx*sx+fy*sy                  # tangential component
    np.fill_diagonal(A, 0.5)              # fill diagonal with 1/2

    # construct the vector `b`
    b = -(np.cos(alpha)*sx+np.sin(alpha)*sy)
```

That seems like the easy part. After all, this is a dense linear algebra problem. How in the world will we determine `gamma`? 

In fact, this is trivial:
```python
    gamma = np.linalg.solve(A, b)
```
##### Numerical fundamental: Linear Algebra Packages
##### Never write your own matrix solver

Every worthwhile numerical language has a set of linear algebra solution routines like the [linalg package](http://docs.scipy.org/doc/numpy/reference/routines.linalg.html). Use them!

---

The function `PanelArray.solve_gamma(alpha=0)` combines the construction and solution code above. After defining __any__ `PanelArray` and an angle of attack `alpha`, this function sets $\gamma_i$ on each panel such that the no-slip condition is enforced.

Let's put all the `VortexPanel` functions together and test them out!


```python
circle = vp.make_circle(N=10)   #1) define geometry
circle.solve_gamma()            #2) solve for gamma
circle.plot_flow()              #3) compute flow field and plot
```

Much better! The flow is generally going around the shape. And the flow is (pretty much) stationary within the body. But the external flow doesn't look __exactly__ like the potential flow solution...

##### Numerical Fundamental: Convergence with resolution
##### The numerical solution is improved by using more panels

---

## Quantitative testing

Let's make this explicitly clear by plotting the distance around the shape $s$ versus $\gamma$ for a number of resolutions.

The function `PanelArray.get_array` lets you get any attribute from the panels. What attributes are there?


```python
help(vp.Panel)
```

So each panel knows all about its geometry (size, location, orrientation) as well as the strength $\gamma$.

Let's try getting this information. Say, the half-width of each panel:


```python
print(circle.get_array('S'))
```

The `distance` function uses this to get the cumulative distance.


```python
print(getsource(circle.distance))
```

The two functions `get_array` and `distance` are all we need to make a quantitative plot of our results:


```python
# Exact solution
s = np.linspace(0,2*np.pi)
plt.plot(s,2*np.sin(s),'k',label='exact')

# Loop over resolutions
for N in 2**np.arange(6,2,-1):         # N in powers of 2
    circle = vp.make_circle(N)         # define geometry
    s = circle.distance()              # get distance array

    circle.solve_gamma()                # solve for gamma
    gamma = circle.get_array('gamma')   # get gamma array
    plt.plot(s,gamma,'--',label=N)      # plot

# finish gamma(s) plot
plt.legend(title='N')
plt.xlabel(r'$s$', fontsize=20)
plt.ylabel(r'$\gamma$', fontsize=20, rotation=0)
plt.show()
```

All of the results for $N>32$ panels look pretty much identical, meaning our results __converge__ as we increase the number of panels. But we can't judge if the converged solution is __correct__ unless we compare to the known result.

##### Numerical fundamental: Validation
##### A method is validated by comparing the result to a known solution

That's easy in this case. The exact solution from analytic potential flow around a circle is

$$\tilde\gamma = -\tilde u_e = 2\sin\theta $$

which is also in the plot and matches the converged solution perfectly.

## Other shapes

Notice that there was nothing special to the circle in this set-up. By giving a different set of points to `panelize`, we can make any shape we need.

##### Quiz

This vortex panel method can be used to solve for the flow around:

1. an ellipse
1. a pair of tandem bodies
1. a rudder

---
##### Your turn

 - ** Modify ** the `make_ellipse` function below to generate an ellipse instead of a circle when supplied with an aspect ratio `t_c`=$t/c$.
 - ** Discuss ** whether the maximum speed around the ellipse is greater or less than that around the circle.
 - ** Move ** a 2:1 ellipse geometry by giving it a different center
 - ** Combine ** the circle and ellipse geometry together into one set of panels using `vp.concatenate` and solve for the flow.
 - ** Discuss ** if there is a *wake* between the bodies. Why or why not?

---

 
##### Solution


```python
def make_ellipse(N, t_c, xcen=0, ycen=0):
    theta = np.linspace(0, -2*np.pi, N+1)
    # your code here to define the points for an ellipse
    x = xcen+np.cos(theta)  # adjust?
    y = ycen+np.sin(theta)  # adjust?
    return vp.panelize(x,y)

# ellipse = make_ellipse(N=32,t_c=2,xcen=2)     # make the shape
# ellipse.solve_gamma()                         # solve for gamma
# ellipse.plot_flow(size=4)                     # compute flow field and plot
```


```python
# pair = ?   your code using vp.concatenate(a1,a2)
```

d174fb6021cb18124b48ed698bc1fa9dc6d8d29c
ipynb
Jupyter Notebook
pre2021/notebooks/3_1_SolutionMethod.ipynb
Maselko/MarineHydro-Tesing
7620a492b629ffad59b30ab300cfb4032f125be5
CC-BY-4.0
2015-03-13T18:42:18.000Z
2021-11-30T11:30:11.000Z
pre2021/notebooks/3_1_SolutionMethod.ipynb
Maselko/MarineHydro-Tesing
7620a492b629ffad59b30ab300cfb4032f125be5
CC-BY-4.0
2020-05-06T12:57:35.000Z
2022-02-15T16:28:55.000Z
pre2021/notebooks/3_1_SolutionMethod.ipynb
Maselko/MarineHydro-Tesing
7620a492b629ffad59b30ab300cfb4032f125be5
CC-BY-4.0
2015-04-10T11:40:13.000Z
2022-01-22T12:42:03.000Z
Qwen/Qwen-72B
1. YES
2. YES