Chapter 3: Potential Flow Theory

mi508668

3 Chapter 3: Potential Flow Theory

Potential flow refers to the movement of a fluid (such as water or air) that relies on assumptions that are consistent with no viscosity or turbulence. This type of flow is often used as an idealized model to describe the behavior of fluids in certain situations, such as in fluid dynamics and aerodynamics. There are incompressible and compressible forms of the model.

In potential flow, the fluid moves in streamlines, which are lines that are tangent to the velocity vector at any given point. The fluid’s velocity is determined by the gradient of the flow’s scalar potential, the flow is irrotational, and the fluid’s pressure is yielded from the Bernoulli equation.

A key concept of potential flow is the construction of a complex flow scenario created by combining different flow models. Some examples of these flow models include:

Source and sink flow: This occurs when at a point (or line in 2D) fluid is injected (source) or ingested (sink) from the flow. The resulting fluid flow will respond. Note that in aerodynamics we often use a “source sheet,” which is a continuous distribution of these. As we will see in class, these are often used to drive models of body thickness.
Doublet flow: This occurs when a source and sink are placed virtually near each other. If the source and sink are coincident, they cancel. However, when very near they create a path of streamlines that form another scenario that looks like a cylinder (2D) and sphere (3D).
Vortex flow: This occurs when there is a spinning motion in the fluid, such as a tornado or a whirlpool without the viscous core. The fluid will flow in a circular pattern around the center of rotation.
Uniform flow: This occurs when the fluid is flowing at a constant velocity in a straight line.

It’s important to note that this type of flow is an idealization and only holds in situations where there is no viscosity and no turbulence in the fluid, and it’s mostly used as a theoretical underpinning for aerodynamic models.

Potential Flow Model

Prior to developing the potential flow model, we must rework how to describe velocity vectors. Potential flow demands defining two scalar fields, the velocity potential function ( $\phi$ ) and the stream function ( $\psi$ ).

Velocity potential, $\phi$ , is 2D and 3D concept, and is, conceptually, a scalar that directly relates to velocity through the gradient of $\phi$ . That is, velocity always moves in the direction of increasing $\phi$ . In this context we obtain the Cartesian form as

(1) $\begin{equation*} \overline{V}=\frac{\partial\phi}{\partial x}\overline{i}+\frac{\partial\phi}{\partial y}\overline{j}+\frac{\partial\phi}{\partial z}\overline{k} \end{equation*}$

Which can be described in vector form as $\overline{V}=\nabla\phi$ . A simple analogy for $\phi$ may be elevation on a topographic map representing a mountain or a hill. The elevation is a scalar. However, if you have a ball it will fall down the steepest path (i.e., the lowest elevation gradient, so the analogy is similar but not exact).

Alternatively, the Stream function, $\psi$ is valid only in 2D, so it only applies to 2D bodies such as airfoils, cylinders, etc. For $\psi$ , a constant value yields a streamline (or a line a particle follows with slopes consistent with the local velocity). We obtain $\psi$ through tangency constraints, or

(2) $\begin{equation*} u=\frac{\partial\psi}{\partial y} \end{equation*}$

(3) $\begin{equation*} v=-\frac{\partial\psi}{\partial x} \end{equation*}$

Now with these definitions, we can move into direct correspondence with the conservation equations.
First, consider the conservation of mass and how to get to potential flow. In the limit of a steady incompressible flow, mass conservation is as follows:

(4) $\begin{equation*} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0 \end{equation*}$

Via the definition of velocity potential, $\overline{V}=\frac{\partial\phi}{\partial x}\overline{i}+\frac{\partial\phi}{\partial y}\overline{j}+\frac{\partial\phi}{\partial z}\overline{k}$ , we can simply insert the definition and we get

(5) $\begin{equation*} \frac{\partial}{\partial x}\frac{\partial\phi}{\partial x}+\frac{\partial}{\partial y}\frac{\partial\phi}{\partial y}+\frac{\partial}{\partial z}\frac{\partial\phi}{\partial z}=0=\nabla^2\phi. \end{equation*}$

Where the Laplacian is defined as follows

(6) $\begin{equation*} \nabla^2\phi=\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}+\frac{\partial^2\phi}{\partial z^2} \end{equation*}$

and conservation of mass in our potential flow model can be represented as

(7) $\begin{equation*} \nabla^2\phi=0. \end{equation*}$

This is known as a Laplace equation which has some of the simplest solutions for a partial differential equation.

In polar coordinates, Laplace equation is given as:

(8) $\begin{equation*} \nabla^2\phi=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial \phi}{\partial r}\right)+\frac{1}{r^2}\frac{\partial^2\phi}{\partial \theta^2}=0 \end{equation*}$

which is important for various flows of interest.

We then bring the Conservation of Momentum back into the context of our Potential Flow model as a scalar relation in the form of the Bernoulli Equation:

(9) $\begin{equation*} \frac{\partial\phi}{\partial t}+\frac{p}{\rho}+\frac{1}{2}V^2+gz=const. \end{equation*}$

Recall that in the context of potential flow, Bernoulli Equation can extend beyond streamlines. Hence, it develops a path to recover pressure if needed for pressure integration.
In the context of potential flow, any equation that satisfies Laplace’s equation is virtually a valid solution for a flow in the limit of potential flow assumptions. The remaining gap is relating this to complex physical scenarios that are relevant to aerodynamic objects of interest, i.e., airfoils, wings, fuselages, etc.
The other aspect of potential flow is the irrotational flow condition, which enables a similar methodology in the context of the stream function. Recall that Irrotational flow is established by the state that the curl is equal to $0$ , or

(10) $\begin{equation*} \nabla\times\overline{V}=0. \end{equation*}$

In the context of a 2D flow, we can expand this to

(11) $\begin{equation*} \nabla\times\overline{V}=\overline{k}\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)=0. \end{equation*}$

Focusing only on the $\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ term and Incorportating the defintition of $u$ and $v$ in terms of the stream function yields

(12) $\begin{equation*} \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}=\frac{\partial}{\partial x}\left(\frac{\partial\psi}{\partial x}\right)-\frac{\partial}{\partial y}\left(-\frac{\partial\psi}{\partial y}\right)=\frac{\partial^2\psi}{\partial x^2}+\frac{\partial^2\psi}{\partial y^2}=0. \end{equation*}$

Hence, the irrotational condition (fundamental to potential flow) demands that the Laplacian of $\psi$ must also be satisfied.

Elementary flows

Elementary potential flows refer to irrotational, inviscid and incompressible fluid flows that can are generated using simple mathematical functions (or flow elements). Examples of elementary potential flows include uniform flow (free stream), source/sink flow, vortex, and doublet. These flows can be used to model the behavior of real-world flows in certain conditions and provide a useful foundation for more complex fluid dynamics analysis. A very nice summary of these is provided in these tables of elementary flows which you should all have on hand. It is important to note these can be written in Cartesian and polar coordinates. As we will find in class, polar is useful for many scenarios.

Superposition

A critical concept for the path to developing complex aerodynamic configurations involves the concept of superimposition of individual potential flow solutions that are simple (elementary or flow elements) to construct models of more complex scenarios. There is a similar analogy with Legos. Each building block has various shapes that individually appear as blocks, but when combined they can recreate more elaborate shapes. In potential flow, our “Lego” blocks are fundamental solutions for (1) free stream flow, (2) source/sink, and/or (3)line vortex (amongst others). We add either a few, or many of these, together to formulate complex models that satisfy the governing equations.

Learning Objectives: Visual description of potential flow

Examples: Superposition

Question: Why can we simply add these simple models together and not violate the governing equation?

If the solution of an individual elementary flow (free-stream, vortex, source) satisfies Laplace’s equation, the result must be zero. It is a trivial solution! Since 0+0=0, we can easily see why superposition works.

Let’s test this out by adding a sink with a vortex at the origin:

(13) $\begin{equation*} \phi=\phi_{Sink}+\phi_{Vortex}=\frac{-m}{2\pi L}\ln{\left(r\right)}+\frac{\Gamma}{2\pi}\theta \end{equation*}$

To evaluate mass conservation (in terms of velocity potential for polar coordinates, Eq. 7):

(14) $\begin{equation*} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\phi}{\partial r}\right)+\frac{1}{r^2}\frac{\partial^2\phi}{\partial\theta^2}=0 \end{equation*}$

Which becomes

(15) $\begin{equation*} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\left(\frac{-m}{2\pi L}\ln{\left(r\right)}+\frac{\Gamma}{2\pi}\theta\right)\right)+\frac{1}{r^2}\frac{\partial^2}{\partial\theta^2}\left(\frac{-m}{2\pi L}\ln{\left(r\right)}+\frac{\Gamma}{2\pi}\theta\right)=0 \end{equation*}$

Yielding

(16) $\begin{equation*} \left(0\right)+\left(0\right)=0. \end{equation*}$

Which works! This can be done for any of these functions. Try on your own for:

A Free Stream
A Free Stream + Vortex @ Origin
A Free Stream + Vortex @ Origin+ Dipole (doublet) @ Origin

An example of this is given in potentialflow.com.

In the context of $\phi$ and $\psi$ it should be clear that the superposition principles apply (e.g., we simply add $\phi$ of various elementary flows). It is also important to note that superposition similarly applies to velocities. That is, the addition of $\phi$ and $\psi$ is equivalent to the addition of the corresponding velocity fields ( $\overbar{V}$ ). It is important to note that $\phi$ and $\psi$ are trivial scalar additions, while $\overbar{V}$ involves vector addition. To explore why consider an example. Let’s assume we have elements 1 and 2 with velocity and velocity potential values given as

(17) $\begin{equation*} u_1=\frac{\partial\phi_1}{\partial x} \end{equation*}$

(18) $\begin{equation*} u_2=\frac{\partial\phi_2}{\partial x} \end{equation*}$

We can obtain $\phi$ as

(19) $\begin{equation*} \phi=\int\ udx \end{equation*}$

Or (note that $u=\frac{\partial\phi}{\partial x}$ , hence, $\phi=\int{udx}+C$ )

(20) $\begin{equation*} \phi_{sum}=\phi_1+\phi_2=\int\ u_1dx+\int\ u_2dx=\int{\left(u_1+u_2\right)dx} \end{equation*}$

It should be observed that the mathematics of addition within the integrals (linking velocity to $\psi$ ) is virtually the same as working with addition in an integrand. Therefore, the addition of $\phi$ (and $\psi$ ) is virtually the same as adding the velocity from each respective elementary flow.

Key Takeaways

Laplace’s equation of $\phi$ is the conservation of mass in the limit of incompressible flow
Laplace’s equation of $\psi$ implies an irrotational flow
Bernoulli equation is one path to understanding pressure (and ultimately integrated to obtain force)
Superposition enables combining the elementary flows, to form more complex flows.
Superposition can be obtained from combining $\phi$ , $\psi$ , or $\overbar{V}$ .

Potential Flow Construction

In order to understand potential flow models, let’s explore various parameters and how they related to altering velocity in a given flow field. For now, lets solely focus on a uniform flow and the addition of sources.

Impact of a Source

Consider a simple scenario: Free Stream + source at origin $(0,0)$ . How does the velocity vector change at $(0,5)$ ?

As can be observed in the figure below, the “free stream” has a velocity ( ${\overline{V}}_\infty$ ) and velocity potential ( $\phi_{FS}$ ) that can be added everywhere in space to that source velocity ( ${\overline{V}}_{Source}$ ) and velocity potential ( $\phi_{Source}$ ). Note that the free-stream velocity, is independent of space (but $\phi_{FS}(x)$ ). For the source, they tend to be a function of $\overline{r}$ from the origin of the source, or ${\overline{V}}_{Source}\left(r\right)$ and $\phi_{Source}(r)$ . For convenience, at $x=0$ , $y=r$ to enable simple first examples. If we focus on the resulting velocity vector, we obtain a total vector that is the resultant vector of the free-stream and source components (as indicated in the figure). Note that the addition of $\phi$ (and $\psi$ ) are scalar additions, hence, do not demand vector addition. What you should focus on is how the total vector can change as we vary the parameters. If we increase the speed, $U_\infty$ , the total vector leans more towards $\parallel$ to the x-axis (and slowing down does the opposite). Alternatively, we can change the source strength, m, to also change this vector orientation. Lastly, although we are fixed in space here, it is important to note that r is also important to change this velocity vector. Let us dig into the relevant details more to better understand.

If you are not following the concept, maybe a first cut could include exploring this flow in the potential flow simulator.

Alternatively, let’s explore analytically. Note for this case, the radial velocity from the source adds a vertical (or $j$ -direction) component. This simplifies the process.

Velocity Potential: Velocity potential ( $\phi$ ) is a scalar, and its addition is trivial. At $(0,5)$ we simply sum the velocity potentials at $(0,5)$ to obtain the combined value:

(21) $\begin{equation*} \phi\left(0,5\right)=U_\infty\left(x\right)+\frac{m}{2\pi}lny\ \end{equation*}$

To obtain velocity we could calculate:

(22) $\begin{equation*} \overline{V}\left(0,5\right)=\left[\frac{\partial\phi}{\partial\ x}\overline{i}+\frac{\partial\phi}{\partial\ y}\overline{j}\right]_{@(0,5)}\ \end{equation*}$

In general, it can be more complex because the source is defined radially, but here

(23) $\begin{equation*} \overline{V}\left(0,5\right)=\left[U_\infty\overline{i}+\frac{m}{2\pi\ y}\ \overline{j}\right]_{@(0,5)} \end{equation*}$

As we will see, this is identical to adding velocity vectors.
Velocity Vector: Velocity vectors are also additive, and their addition is not too difficult, especially in this scenario. At $(0,5)$ we again sum the velocity vector contributions at $(0,5)$ to obtain the combined value. To obtain velocity we could calculate:

(24) $\begin{equation*} \overline{V}(0,5)=U_\infty\overline{i}+\frac{m}{2\pi(5m)}\overline{j}\ \end{equation*}$

However, this is complex because the source is defined radially. So the velocity vector becomes a function of m (which has units $m^2/s$ ). For a source, $m\gt 0$ , and the velocity points upward. Alternatively, for a sink m<0 and the velocity will point downward. As an example, if we want a vector to have a $45 deg$ angle with respect to the x-axis, we need to match the $i$ and $j$ components, or:

(25) $\begin{equation*} U_\infty=\frac{m}{2\pi(5m)}\rightarrow\ m=2\pi\left(5m\right)\left(\frac{100m}{s}\right) \end{equation*}$

This is used to establish flow tangency with walls in a potential flow model.

Examples: Sample Matlab Script

clear

pi=2*acos(0);

%uniform Flow

Uinf=100;

%Source

xS=0;

yS=0;

%Point of interest

xP=0;

yP=5.0;

%solve for source strengh that provides a 45 deg angle (u=v@(0,5))

uP=Uinf;

vP=Uinf;

m=vP*(2*pi*sqrt((yP-yS)^2+(xP-xS)^2));

x=-10:10;y=-10:10;

[X,Y]=meshgrid(x,y)

R=sqrt((X-xS).^2+(Y-yS).^2);%compute R

THETA=atan2(Y,X);%compute THETA

%build phi

phi=Uinf*X+m/(2*pi)*log(R);

%build psi

psi=Uinf*Y+m/(2*pi)*THETA;

%build Vs

dbar(:,:,1)=(X-xS)./R;

dbar(:,:,2)=(Y-yS)./R;

%find Velocity

Vel(:,:,1)=Uinf+dbar(:,:,1)*m./(2*pi*abs(R));

Vel(:,:,2)=dbar(:,:,2)*m./(2*pi*abs(R));

%plot results

contourf(X,Y,phi)

hold on

contour(X,Y,psi)

quiver(X,Y,Vel(:,:,1),Vel(:,:,2))

hold off

Key Takeaways

In this simple example, we found that:

The free stream acts like a free stream velocity in aerodynamics. It is a boundary condition.
The source strength and distance control velocity vectors at specific points. This is important as it is used to specify “flow tangency” on surfaces.

Impact of Source Location

Consider a slightly more complex scenario: Free Stream + source at $(5,0)$ . How does the velocity vector change at $(0,5)$ ?

You can try this in the potential flow simulator.

Alternatively, let’s explore analytically. For this case, the radial velocity from the source DOES NOT add a vertical (or $j$ -direction) component to the velocity vector. This can be seen in the vector addition in the figure above. This makes the process a little more complex. But let’s again explore using the velocity potential and vector addition.

Velocity Potential: Again, velocity potential ( $\phi$ ) is a scalar, and its addition is trivial. At $(0,5)$ we simply sum the velocity potentials at $(0,5)$ to obtain the combined value:

(26) $\begin{equation*} \phi\left(0,5\right)=U_\infty\left(x\right)+\frac{m}{2\pi}lnr=\frac{m}{2\pi}\ln{\left(\sqrt{50}m\right)}\ \ \end{equation*}$

This scalar quantity is easy to work with but we have the same issue with obtaining velocity (mixing Cartesian and radial coordinates)

Velocity Vector: Velocity vectors are still additive, but demand to be accounted for in mixed coordinate systems. Similar to before, at $(0,5)$ we again sum the velocity vector contributions at $(0,5)$ to obtain the combined value. This time, however, the source acts in both the $i$ and $j$ directions. The magnitude is given by the radius

(27) $\begin{equation*} {|V}_s(0,5)|=\frac{m}{2\pi\left|r\right|} \end{equation*}$

The radial direction considered here demands decomposing into $i$ and $j$ components. This demands defining the vector $\overline{r}$ as indicated in the figure. This vector provides direction in Cartesian coordinates. The direction is given by:

(28) $\begin{equation*} \hat{d}=\frac{\overline{r}}{\left|r\right|}=-\frac{5}{\sqrt{50}}\overline{i}+\frac{5}{\sqrt{50}}\overline{j} \end{equation*}$

Just a note, when we compute the velocity from a source at point $p$ (which is $(0,5)$ in this example), we can generalize as follows:
$\overline{r}={\overline{r}}_p=\left(x_p-x_{source}\right)\overline{i}+(y_p-y_{source})\overline{j}$
We can now sum the velocity components of the free stream and source as follows:

(29) $\begin{equation*} \overline{V}\left(0,5\right)=U_\infty\overline{i}+\frac{m}{2\pi\left|r\right|}\frac{\overline{r}}{\left|r\right|}=\left(100+\frac{m}{2\pi\sqrt{50}}\frac{\left(-5\right)}{\sqrt{50}}\right)\overline{i}+\ \left(\frac{m}{2\pi\sqrt{50}}\frac{5}{\sqrt{50}}\right)\overline{j} \end{equation*}$

As an example, if we want a vector to have a $45$ deg angle with respect to the x-axis, we need to match the $i$ and $j$ components, or:

(30) $\begin{equation*} \left(100-\frac{m}{2\pi\sqrt{50}}\frac{5}{\sqrt{50}}\right)=\left(\frac{m}{2\pi\sqrt{50}}\frac{5}{\sqrt{50}}\right)\rightarrow\left(100\right)=\frac{m}{2\pi\sqrt{50}}\left(\frac{10}{\sqrt{50}}\right) \end{equation*}$

The result yields $m=1000\pi$ . Hence, flow tangency with walls in a potential flow model can be constructed from sources in more complex configurations.

Examples: Matlab Script of uniform flow and source @ (5,0)

clear

pi=2*acos(0);

%uniform Flow

Uinf=100;

%Source

xS=5;

yS=0;

%Point of interest

xP=0;

yP=5.0;

%solve for source strengh that provides a 45 deg angle (u=v@(0,5))

uP=Uinf;

vP=Uinf;

m=vP*(2*pi*sqrt((yP-yS)^2+(xP-xS)^2));

x=-10:10;y=-10:10;

[X,Y]=meshgrid(x,y)

R=sqrt((X-xS).^2+(Y-yS).^2);%compute R

THETA=atan2(Y,X);%compute THETA

%build phi

phi=Uinf*X+m/(2*pi)*log(R);

%build psi

psi=Uinf*Y+m/(2*pi)*THETA;

%build Vs

dbar(:,:,1)=(X-xS)./R;

dbar(:,:,2)=(Y-yS)./R;

%find Velocity

Vel(:,:,1)=Uinf+dbar(:,:,1)*m./(2*pi*abs(R));

Vel(:,:,2)=dbar(:,:,2)*m./(2*pi*abs(R));

%plot results

contourf(X,Y,phi)

hold on

contour(X,Y,psi)

quiver(X,Y,Vel(:,:,1),Vel(:,:,2))

hold off

Key Takeaways

A source induces a velocity in the radial direction which must be considered from the velocity “induced” by the source.

Multiple Sources

Consider a more complex scenario: Free Stream + source at $(0,0)$ + source at $(5,0)$ . How does the velocity vector change at $(0,5)$ ? The scenario is indicated in the Figure below and can be replicated in the potential flow simulator.

Again, we need to explore this scenario analytically. For this case, we need to add components from multiple sources. This can be seen in the vector addition in the figure above. This process is leading to how a panel code works. Lets us again explore using the velocity potential and vector addition.

Velocity Potential: Again, velocity potential ( $\phi$ ) is a scalar, and its addition is trivial. At $(0,5)$ we simply sum the velocity potentials at $(0,5)$ to obtain the combined value:

(31) $\begin{equation*} \phi\left(0,5\right)=U_\infty\left(x\right)+\frac{m_1}{2\pi}lnr_1+\frac{m_2}{2\pi}lnr_2=\frac{m_1}{2\pi}\ln{\left(5m\right)}+\frac{m_2}{2\pi}\ln{\left(\sqrt{50}m\right)}\ \ \end{equation*}$

We continue to see the ease of working with this scalar quantity. But we also need to work with velocities.
Velocity Vector: Velocity vectors remain additive and this time include mixed coordinate systems and multiple sources. Similar to before, at $(0,5)$ we again sum the velocity vector contributions at $(0,5)$ to obtain the combined value. This time, however, two sources acts and need to be additively evaluated. When we sum the velocity components of the free stream and two source flows, we obtain:

(32) $\begin{equation*} \overline{V}\left(0,5\right)=U_\infty\overline{i}+\frac{m_1}{2\pi\ r_1}\frac{\overline{r_1}}{\left|r_1\right|}+\frac{m_2}{2\pi\ r_2}\frac{\overline{r_2}}{\left|r_2\right|} \end{equation*}$

which becomes

(33) $\begin{equation*} \overline{V}\left(0,5\right)=\left(100+\frac{m_1}{2\pi\ r_1}\frac{(0)}{\sqrt5}-\frac{m_2}{2\pi\ r_2}\frac{5}{\sqrt{50}}\right)\overline{i}+\ \left(\frac{m_1}{2\pi\ r_1}\frac{5}{\sqrt5}+\frac{m_2}{2\pi\ r_2}\frac{5}{\sqrt{50}}\right)\overline{j} \end{equation*}$

This can be expanded into a matrix multiplication as follows:

(34) $\begin{equation*} \overline{V}\left(x,y\right)={\overline{A}}^T\overline{m}+{\overline{V}}_\infty \end{equation*}$

Where ${\overline{A}}^T$ is a vector of influence coefficients, given as:

(35) $\begin{equation*} {\overline{A}}^T=\frac{1}{2\pi}\left[\begin{matrix}\frac{\overline{r_1}}{\left|r_1\right|^2}&\frac{\overline{r_2}}{\left|r_2\right|^2}\\\end{matrix}\right] \end{equation*}$

There is a vector of source strengths given as:

(36) $\begin{equation*} \overline{m}=\left[\begin{matrix}m_1\\m_2\\\end{matrix}\right] \end{equation*}$

This can be generalized to the following form where there are N_{sources}:

(37) $\begin{equation*} \overline{V}\left(x,y\right)={\overline{V}}_\infty+\sum_{i=1}^{N_{sources}}\frac{1}{2\pi}\frac{\overline{r_i}}{\left|r_i\right|^2}m_i={\overline{V}}_\infty+{\overline{A}}^T\overline{m} \end{equation*}$

The previous examples indicated how to find a source strength to specify a velocity direction. In this case, we have 1 equation and 2 unknowns ( $m1$ and $m2$ ), hence, it is not solvable unless we fix one of the $m's$ . Because of this we need to add a new equation to specify a velocity at another point. Here is a sample matlab/octave script:

Examples Matlabe script for multiple sources

clear

pi=2*acos(0);

%uniform Flow

Uinf=100;

%Source

xS(1)=0;%source 1

yS(1)=0;

xS(2)=5;%source 2

yS(2)=0;

%This time just set the source strengths

m(1)=1000;

m(2)=-500

x=-10:10;y=-10:10;

[X,Y]=meshgrid(x,y)

phi=X*0;%initialize phi to all 0;

psi=phi;%initialize psi to all 0;

Vel(:,:,1)=phi; %initialize vel to 0

Vel(:,:,2)=phi;

for ii=1:2 %loop over all sources

%Note that this is all additive and phi,psi, and Vel all build each

%source

R=sqrt((X-xS(ii)).^2+(Y-yS(ii)).^2);%compute R1

THETA=atan2(Y-yS(ii),X-xS(ii));%compute THETA

%build phi

phi=m(ii)/(2*pi)*log(R)+phi;

%build psi

psi=m(ii)/(2*pi)*THETA+psi;

%build Vs

dbar(:,:,1)=(X-xS(ii))./R;

dbar(:,:,2)=(Y-yS(ii))./R;

%find Velocity

Vel(:,:,1)=dbar(:,:,1)*m(ii)./(2*pi*abs(R))+Vel(:,:,1);

Vel(:,:,2)=dbar(:,:,2)*m(ii)./(2*pi*abs(R))+Vel(:,:,2);

end

%add free stream component

phi=Uinf*X+phi;

psi=Uinf*Y+psi;

Vel(:,:,1)=Vel(:,:,1)+Uinf;

%plot results

contourf(X,Y,phi)

hold on

contour(X,Y,psi)

quiver(X,Y,Vel(:,:,1),Vel(:,:,2))

hold off

Key Takeaways

Multiple sources induce additive velocity, each with their own radial direction.
Closing out the source strengths demands 1 equation for each source. This implies you will need to enforce velocity at two locations to obtain two sources.

Constructing a “no penetration” or “flow tangency” condition

Consider a more complex scenario: Free Stream + source at (0,5) + source at (5,0) and flat plate walls angled of +45 degrees (0,0) and -45 deg (5,5). The scenario is depicted in the figure below. Although the concept is complex, the underlying mathematics is not overly onerous.

Figure: two source/two flow tangency scenario.

Now lets move into an approach to find source strengths. Recall from last time, we had two unknowns (m1 and m2) and one equation (velocity at (0,5)). This time, we have two flow tangency conditions giving us two equations (to close out the two unknowns). This time let’s just focus on the velocity vectors. Evaluating on wall 1 (w1) we get:

(38) $\begin{equation*} \overline{V}\left(0,0\right)={\overline{V}}_{w1}={\overline{V}}_\infty+\sum_{j=1}^{N_{sources}}\frac{1}{2\pi}\frac{\overline{r_{w1j}}}{\left|r_{w1j}\right|^2}m_i \end{equation*}$

which yields

(39) $\begin{equation*} \overline{V}\left(0,0\right)={\overline{V}}_\infty+\frac{1}{2\pi}\frac{\overline{r_{w11}}}{\left|r_{w11}\right|^2}m_1+\frac{1}{2\pi}\frac{\overline{r_{w12}}}{\left|r_{w12}\right|^2}m_2. \end{equation*}$

Note the convention, $r_{ij}$ is the radius from a source at location $j$ acting at point (or wall location) $i$ . This equation has two unknowns, $m1$ and $m2$ . At $w1$ we seek to find flow tangency, that is flow is parallel to the wall. This is invoked in the condition where:

(40) $\begin{equation*} {\overline{V}}_{w1}\cdot{\overline{n}}_1=0 \end{equation*}$

Expanding this we obtain:

(41) $\begin{equation*} {\overline{V}}_{w1}\cdot{\overline{n}}_1=\left({\overline{V}}_\infty+\frac{1}{2\pi}\frac{\overline{r_{11}}}{\left|r_{11}\right|^2}m_1+\frac{1}{2\pi}\frac{\overline{r_{12}}}{\left|r_{12}\right|^2}m_2\right)\cdot{\overline{n}}_1=0 \end{equation*}$

Or

(42) $\begin{equation*} {\overline{V}}_\infty\cdot{\overline{n}}_1+\frac{1}{2\pi}\frac{\overline{r_{11}}\cdot{\overline{n}}_1}{\left|r_{11}\right|^2}m_1+\frac{1}{2\pi}\frac{\overline{r_{12}}\cdot{\overline{n}}_1}{\left|r_{12}\right|^2}m_2=0 \end{equation*}$

Performing the same procedure on w2, we obtain:

(43) $\begin{equation*} {\overline{V}}_\infty\cdot{\overline{n}}_2+\frac{1}{2\pi}\frac{\overline{r_{21}}\cdot{\overline{n}}_2}{\left|r_{21}\right|^2}m_1+\frac{1}{2\pi}\frac{\overline{r_{22}}\cdot{\overline{n}}_2}{\left|r_{22}\right|^2}m_2=0 \end{equation*}$

At this point, we have two equations and 2 unknowns, hence a solvable (and linear) system is formed. For this smaller case, a manual solution is possible. However, what we really seek it to extend to more complex and larger systems. Hence, we introduce a matrix form. We can combine the proceeding two equations into matrix in terms of $\left[A\right]\overline{x}=\overline{b}$ given as

(44) $\begin{equation*} \left[A\right]=\frac{1}{2\pi}\begin{bmatrix} \frac{\bar{r}_{11}\cdot\bar{n}_{1}}{\left|\bar{r}_{11}^2 \right|} & \frac{\bar{r}_{12}\cdot\bar{n}_{1}}{\left|\bar{r}_{12}^2 \right|} \\ \frac{\bar{r}_{21}\cdot\bar{n}_{2}}{\left|\bar{r}_{21}^2 \right|} & \frac{\bar{r}_{22}\cdot\bar{n}_{2}}{\left|\bar{r}_{22}^2 \right|} \end{bmatrix}. \end{equation*}$

Here the elements of $\left[A\right]$ are called influence coeffs and given as

(45) $\begin{equation*} A_{i,j}=\frac{\overline{r_{ij}}\cdot{\overline{n}}_i}{\left|r_{ij}\right|^2}\ \end{equation*}$

where the elements are given as $x_i=m_i$ and

(46) $\begin{equation*} \bar{b}=-\left[\begin{matrix}{\bar{V}}_\infty\cdot{\overline{n}}_1\\{\bar{V}}_\infty\cdot{\overline{n}}_2\\\end{matrix}\right] \end{equation*}$

where the elements are given as $b_i=-{\overline{V}}_\infty\cdot{\overline{n}}_i$ . It is important that you see how this system is formed. Please review and make sure you understand it.

We can now solve using linear algebra techniques. Fortunately, we do not need fancy ones, and a brute force, matrix inversion, will do the trick. That is

(47) $\begin{equation*} \overline{x}=\left[A\right]^{-1}\overline{b} \end{equation*}$

In this solution, $\overline{x}$ is a vector of source strengths that will simultaneously satisfy flow tangency at the two constraints.

Examples: Sample matlab code to determine source strength to satisfy flow tangency

clear

pi=2*acos(0);

%uniform Flow

Uinf=100;

%Source

xS(1)=0;%source 1

yS(1)=5;

xS(2)=5;%source 2

yS(2)=0;

%walls

xW(1)=0; %location

yW(1)=0;

nw(1,1)=-1/sqrt(2); %wall 1, i dir

nw(1,2)= 1/sqrt(2); %wall 1, j dir

xW(2)=5; %location

yW(2)=5;

nw(2,1)= 1/sqrt(2); %wall 2, i dir

nw(2,2)= 1/sqrt(2); %wall 2, j dir

%This time we need to find the source strengths

%create A matrix (that is Aij)

% Note i is the row and correspones to the wall

% j is the column and corresponds to the source

for ii=1:2 %loop over each wall

for jj=1:2 %loop over each source

%compute R (so vector form source to point or wall)

rij(1)=xW(ii)-xS(jj);

rij(2)=yW(ii)-yS(jj);

rijMag=sqrt(rij(1)^2+rij(2)^2);

A(ii,jj)=1/(2*pi)*dot(rij,nw(ii,:))/rijMag^2;

end

%we can also load the b vector here

b(ii)=-Uinf*nw(ii,1); %this could be a dot product but Vinf is only in x

end

m=A^-1*transpose(b); %need to transpose b to make a single column vect.

%Setup mesh for flow output

x=-10:10;y=-10:10;

[X,Y]=meshgrid(x,y)

phi=X*0;%initialize phi to all 0;

psi=phi;%initialize psi to all 0;

Vel(:,:,1)=phi; %initialize vel to 0

Vel(:,:,2)=phi;

%This part outputs the flow

for ii=1:2 %loop over all sources

%Note that this is all additive and phi,psi, and Vel all build each

%source

R=sqrt((X-xS(ii)).^2+(Y-yS(ii)).^2);%compute R1

THETA=atan2(Y-yS(ii),X-xS(ii));%compute THETA

%build phi

phi=m(ii)/(2*pi)*log(R)+phi;

%build psi

psi=m(ii)/(2*pi)*THETA+psi;

%build Vs

dbar(:,:,1)=(X-xS(ii))./R;

dbar(:,:,2)=(Y-yS(ii))./R;

%find Velocity

Vel(:,:,1)=dbar(:,:,1)*m(ii)./(2*pi*abs(R))+Vel(:,:,1);

Vel(:,:,2)=dbar(:,:,2)*m(ii)./(2*pi*abs(R))+Vel(:,:,2);

end

%add free stream component

phi=Uinf*X+phi;

psi=Uinf*Y+psi;

Vel(:,:,1)=Vel(:,:,1)+Uinf;

%plot results

contourf(X,Y,phi)

hold on

contour(X,Y,psi)

quiver(X,Y,Vel(:,:,1),Vel(:,:,2))

hold off

Key Takeaways

For a series of unknown sources intended to drive flow tangency conditions. We need to satisfy flow tangency at one point for every source.
The resulting equation set provides a linear system of equations specific to the flow tangency condition.

Source Panels

We can extend these concepts to what is referred to as a source panel. The concepts actually arrive from a basic potential flow concept of a source sheet shown in the diagram below:

Figure: A diagram of a source sheet, which is a continuous distribution of sources that form a “sheet”

What we seek to do is position a “source sheet” on the walls of an aerodynamic body where we also seek to satisfy flow tangency (when also exposed to a free-stream velocity). We also want to simplify and treat the sheet as linear (i.e., line with slope and no curvature). To get to this scenario, consider: Free Stream + 2 m source sheet at (0,5) and positioned at a 0 deg orientation with a width of $1.5 m$ . Then another source at (5,0) and angled at -45 deg (also 1.5m in width). Note that the source sheets are aligned with thin plates. The scenario is depicted below:

In this scenario (depicted below) we need to add two concepts:

Accounting for sheets rather than point sources
Accounting for self-influence (ie, the effect of the source sheet on itself)

To move from a source sheet to a point (or line) source we need to rework our source description. With the sheet we have an integral as follows:

(48) $\begin{equation*} V_r\left(x,y\right)=\int_{s=a}^{b}{\frac{\lambda\left(s\right)}{2\pi\ r(s)}ds}\ \ \ \end{equation*}$

Note that the radius, $r$ , and source/unit length, \lambda, are both a function of the sheet. The goal is to use a first-order method that uses a mean value of $\lambda$ , $\bar{\lambda}$ , evaluated and the center of the source sheet (implying $\bar{r}$ ). These assumptions enable the following simplification.

(49) $\begin{equation*} V_r\left(x,y\right)=\int_{s=a}^{b}{\frac{\lambda\left(s\right)}{2\pi\ r(s)}ds}\cong\frac{\lambda_{avg}}{2\pi\ r_{avg}}\int_{s=a}^{b}ds=\frac{\lambda_{avg}l}{2\pi\ r_{avg}}\ \ \ \end{equation*}$

Such an approximation is provided pictorially below.

Figure: Approximation of source sheet at a single point at the center of the panel.

Using this approach, our velocity equation for an arbitrary point p becomes:

(50) $\begin{equation*} V_r\left(x,y\right)=\int_{s=a}^{b}{\frac{\lambda\left(s\right)}{2\pi\ r(s)}ds}\ \ \ \end{equation*}$

Again, invoking the flow tangency condition:

(51) $\begin{equation*} {\bar{V}}_{w1}\cdot{\bar{n}}_1=0 \end{equation*}$

We obtain the following flow tangency constraints for w1:

(52) $\begin{equation*} {\bar{V}}_\infty\cdot{\bar{n}}_1+\frac{1}{2\pi}\frac{\bar{r_{11}}\cdot{\bar{n}}_1}{\left|r_{11}\right|^2}\lambda_{avg,1}l_1+\frac{1}{2\pi}\frac{\bar{r_{12}}\cdot{\bar{n}}_1}{\left|r_{12}\right|^2}\lambda_{avg,2}l_2=0 \end{equation*}$

and for w2:

(53) $\begin{equation*} {\bar{V}}_\infty\cdot{\bar{n}}_2+\frac{1}{2\pi}\frac{\bar{r_{21}}\cdot{\bar{n}}_2}{\left|r_{21}\right|^2}\lambda_{avg,1}l_1+\frac{1}{2\pi}\frac{\bar{r_{22}}\cdot{\bar{n}}_2}{\left|r_{22}\right|^2}\lambda_{avg,2}l_2=0 \end{equation*}$

In this scenario, we have an issue with $r_{11}$ and $r_{22}$ . The issue is that they both approach zero, leading to a divide by zero. We invoke the

(54) $\begin{equation*} \lim_{r\to0}{\left(\frac{1}{2\pi}\frac{\overline{r}\cdot\overline{n}}{\left|r\right|^2}\lambda_{avg}l\right)\rightarrow\left(\frac{1}{2}\lambda_{avg}\right)\ } \end{equation*}$

Using this approximation of the induced flow on itself, we eliminate the numerical issue. We can combine the proceeding two equations into matrix in terms of $\left[A\right]\bar{x}=\bar{b}$ given as

(55) $\begin{equation*} \left[A\right]=\frac{1}{2\pi}\begin{bmatrix} \pi & \frac{\bar{r}_{12}\cdot\bar{n}_{1}}{\left|\bar{r}_{12}^2 \right|}l_2 \\ \frac{\bar{r}_{21}\cdot\bar{n}_{2}}{\left|\bar{r}_{21}^2 \right|}l_1 & \pi \end{bmatrix}. \end{equation*}$

These influence coefficients are generalized as

(56) $\begin{equation*} A_{i,j}=\frac{\overline{r_{ij}}\cdot{\overline{n}}_il_j}{\left|r_{ij}\right|^2}\ \ for\ i\neq\ j and \pi\ for\ i=j \end{equation*}$

Our source strengths are defined as

(57) $\begin{equation*} \bar{x}=\left[\begin{matrix}\lambda_{avg,1}\\\lambda_{avg,2}\\\end{matrix}\right] \end{equation*}$

where the elements are given as $x_i=\lambda_{avg,i}$ . Finally, the remaining vector does not change and stays as

(58) $\begin{equation*} \bar{b}=-\left[\begin{matrix}{\bar{V}}_\infty\cdot{\bar{n}}_1\\{\bar{V}}_\infty\cdot{\bar{n}}_2\\\end{matrix}\right] \end{equation*}$

where the elements are given as $b_i=-{\bar{V}}_\infty\cdot{\bar{n}}_i$ . Similar to before, we can now solver using linear algebra techniques and solve

(59) $\begin{equation*} $\bar{x}=\left[A\right]^{-1}\overline{b}$. \end{equation*}$

Which yields the final values of the source strengths that satisfies flow tangency.

Key Takeaways

We move from a line source to a source sheet of finite length. This source sheet will become our “panels” that we will later use to represent a surface. We use a mean source strength and radius leading to a 1st order approximation. Anderson presents a higher-order method, but it is not instructive to look at it as a higher order.
We also have to deal with singularity associated with the self-induced flow. The approach to remove this was given. Generally speaking, such singularities challenge potential flow methods

Panel Code

The aforementioned discussion leads to aerodynamic potential flow codes or panel codes. These codes are used to determine velocities and pressure distributions on objects. Note that this is only one path of potential flow. As we will see, it also leads to the foundational theory of aerodynamics. Panel codes can be thought of (at least as a first level) as a series of sources, sinks, vortex points, and/or doublets being positioned along panels that approximately (and discretely) represent an aerodynamic body. These panels can have various lengths and strengths.

There are also equations that can be formulated based on the summation of all the flow contributors (i.e., free stream velocity + induced flow from each panel) and the kinematic constraint of satisfying flow tangency on each panel. This forms a system of linear equations with unknowns that is solvable. This essentially extends the “source panels” section to closed bodies of interest. In this discussion, we aim to develop how sources can be combined with a uniform flow to develop a panel code within a stepwise procedure.

Figure: Sample paneled airfoil at a free-stream velocity with incidence.

Step 1: Define the BCs

The first step demands free stream flight conditions including the speed and angle of attack. We also have to define the shape or paneling of the object. An example of this is depicted in the Figure above. It is important to note that each panel will (by default) be exposed to the free-stream velocity vector. Our goal aims to find the solution that yields flow tangency on each of the panels.

Step 2: Compute the control points

With the paneled airfoil, we now need to compute the “control points” on each panel. This is a point where we can approximately apply the sources and flow tangency. At this level let’s just assume we only need to evaluate at this control point. The control point is given as the midpoint of the panel:

(60) $\begin{equation*} {\bar{x}}_{cp,i}=0.5\left({\bar{x}}_{v,i}+{\bar{x}}_{v,i+1}\right) \end{equation*}$

Note that ${\bar{x}}_{cp,i}$ is the control point (vector) of panel i and ${\bar{x}}_{v,i}$ and ${\bar{x}}_{v,i+1}$ are the vertices of panel i. Note that the panels should be ordered in a consistent direction around the closed body of interest (either clockwise or counterclockwise). This is important to keep the panel number as described correctly. Once defined, we formulate a model that satisfies flow tangency on the control points. The key need is to define the induced velocity from the other panels which involves calculated distances and directions.

Step 3: Compute the panel length and slope or normal vector

Closing this out demands a few other criteria:

Panel length: we solve for an average strength, but it needs to consider being applied to a finite-length source sheet. $ds_i=\left|{\bar{x}}_{v,i+1}-{\bar{x}}_{v,i}\right|$
Panel Slope: Enables us to set up a flow tangency condition (we can also define the normal vector) ${\frac{dy}{dx}}_i=\frac{y_{v,i+1}-y_{v,i}}{x_{v,i+1}-x_{v,i}}$ Or, the (outward) normal vector given as ${\bar{n}}_i=\frac{\left(y_{v,i+1}-y_{v,i}\right)\bar{i}-\left(x_{v,i+1}-x_{v,i}\right)\bar{j}}{\sqrt{\left(x_{v,i+1}-x_{v,i}\right)^2+\left(y_{v,i+1}-y_{v,i}\right)^2}}$

Step 4: Now compute the influence coefficients

This is an extension of the previously discussed equations (in a 2-panel system).

(61) $\begin{equation*} \left[A\right]=\frac{1}{2\pi}\begin{bmatrix} \pi & \frac{\bar{r}_{12}\cdot\bar{n}_{1}}{\left|\bar{r}_{12}^2 \right|}l_2 \\ \frac{\bar{r}_{21}\cdot\bar{n}_{2}}{\left|\bar{r}_{21}^2 \right|}l_1 & \pi \end{bmatrix} \end{equation*}$

where the elements (called influence coefficients) and given as

(62) $\begin{equation*} A_{i,j}=\frac{\bar{r_{ij}}\cdot{\bar{n}}_i}{\left|r_{ij}\right|^2}\ \bar{x}=\left[\begin{matrix}m_1\\m_2\\\end{matrix}\right] \end{equation*}$

where the elements are given as $x_i=m_i$ and

(63) $\begin{equation*} \bar{b}=-\left[\begin{matrix}{\overline{V}}_\infty\cdot{\overline{n}}_1\\{\bar{V}}_\infty\cdot{\bar{n}}_2\\\end{matrix}\right] \end{equation*}$

where the elements are given as $b_i=-{\bar{V}}_\infty\cdot{\bar{n}}_i$ .

Step 5: Reconstruct velocity and pressure coefficient

This is a vector approach that demands creating a function. Then $C_p=1-\left(\frac{\left|V\left(x,y\right)\right|}{|V_\infty|}\right)^2$

Examples: Panel Code Walkthrough

3 Chapter 3: Potential Flow Theory

Potential Flow Model

Elementary flows

Superposition

Potential Flow Construction

Impact of a Source

Impact of Source Location

Multiple Sources

Constructing a “no penetration” or “flow tangency” condition

Source Panels

Panel Code

Step 1: Define the BCs

Step 2: Compute the control points

Step 3: Compute the panel length and slope or normal vector

Step 4: Now compute the influence coefficients

Step 5: Reconstruct velocity and pressure coefficient

License

Share This Book