Version 0.1 -- May 2014¶
2D Source-vortex panel method with boundary-layer correction¶
A source-panel method leads to a solution with no circulation, therefore no lift. The objective of this lesson is to start with the source panel method we implemented in the previous lesson and add some circulation so that we may have a lift force. We introduce an important concept: the Kutta-condition that allows us to determine what the right amount of circulation should be.
In all prior notebooks we have neglected all viscous forces. This notebook will add a boundary layer correction to the potential flow code around a NACA0012 airfoil.
In two dimensional flow around an airfoil one can consider the viscous effects to be confined to a thin layer on the surface. The boundary layer, as coined by Prandtl, is a region just above the surface where the frictional effects, caused by the no-slip condition, effect the flow tremedously. Outside the boundary layer the flow can be considered inviscid, hence potential flow theory applies. The impact of this concept by Prandtl, the boundary layer, on modern aerodynamics cannot be overstated.
The boundary layer begins to grow the moment the flow comes in contact with the surface of the airfoil and continues to do so until the trailing edge. At a certain point in its evolution the boundary layer transitions from a laminar scheme to a turbulent one. To apply these concepts to the code we've been working with we will employ three methods: Thwaites', Michaels, and Heads. Thwaites' method applies to the laminar boundary layer, Michael's criterion determines the transition point, and Head's method applies to the turbulent boundary layer.
Discretizition of airfoil geometry¶
Let's get the preliminaries out of the way. We need to import our favorite libraries, and the function integrate from SciPy, as in Lesson 10.
from math import *
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
# reads of the geometry from a data file
with open ('../resources/naca0012.dat') as file_name:
x, y = np.loadtxt(file_name, dtype=float, delimiter='\t', unpack=True)
# plots the geometry
%matplotlib inline
val_x, val_y = 0.1, 0.2
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()
x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)
y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)
size = 10
plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
plt.grid(True)
plt.xlabel('x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.plot(x, y, color='k', linestyle='-', linewidth=2);
We define a class Panel that will store all information about one panel: start and end points, center point, length, orientation, source strength, tangential velocity and pressure coefficient. We don't save the vortex-sheet strength because all panels will have the same value.As with both Lesson 10 and 11, we define the attributes of a panel with the class Panel.
class Panel:
"""Contains information related to one panel."""
def __init__(self, xa, ya, xb, yb):
"""Creates a panel.
Arguments
---------
xa, ya -- Cartesian coordinates of the first end-point.
xb, yb -- Cartesian coordinates of the second end-point.
"""
self.xa, self.ya = xa, ya
self.xb, self.yb = xb, yb
self.xc, self.yc = (xa+xb)/2, (ya+yb)/2 # control-point (center-point)
self.length = sqrt((xb-xa)**2+(yb-ya)**2) # length of the panel
# orientation of the panel (angle between x-axis and panel's normal)
if xb-xa <= 0.:
self.beta = acos((self.yb-self.ya)/self.length)
elif xb-xa > 0.:
self.beta = pi + acos(-(self.yb-self.ya)/self.length)
# location of the panel
if self.beta <= pi:
self.side = 'extrados'
else:
self.side = 'intrados'
self.sigma = 0. # source strength
self.vt = 0. # tangential velocity
self.cp = 0. # pressure coefficient
Like before, we call the function define_panels() to discretize the airfoil geometry in N panels. The function will return a NumPy array of N objects of the type Panel.
def define_panels(x, y, N=40):
"""Discretizes the geometry into panels using 'cosine' method.
Arguments
---------
x, y -- Cartesian coordinates of the geometry (1D arrays).
N - number of panels (default 40).
Returns
-------
panels -- Numpy array of panels.
"""
R = (x.max()-x.min())/2 # radius of the circle
x_center = (x.max()+x.min())/2 # x-coord of the center
x_circle = x_center + R*np.cos(np.linspace(0, 2*pi, N+1)) # x-coords of the circle
x_ends = np.copy(x_circle) # projection of the x-coords on the surface
y_ends = np.empty_like(x_ends) # initialization of the y-coords Numpy array
x, y = np.append(x, x[0]), np.append(y, y[0]) # extend arrays using np.append
# computes the y-coordinate of end-points
I = 0
for i in range(N):
while I < len(x)-1:
if (x[I] <= x_ends[i] <= x[I+1]) or (x[I+1] <= x_ends[i] <= x[I]):
break
else:
I += 1
a = (y[I+1]-y[I])/(x[I+1]-x[I])
b = y[I+1] - a*x[I+1]
y_ends[i] = a*x_ends[i] + b
y_ends[N] = y_ends[0]
# creates panels
panels = np.empty(N, dtype=object)
for i in xrange(N):
panels[i] = Panel(x_ends[i], y_ends[i], x_ends[i+1], y_ends[i+1])
return panels
Now we can use our new function to define the geometry for the airfoil panels, and then plot the panel nodes on the geometry.
N = 80 # number of panels
panels = define_panels(x, y, N) # discretizes of the geometry into N panels
# plots the geometry and the panels
val_x, val_y = 0.1, 0.2
x_min, x_max = min( panel.xa for panel in panels ), max( panel.xa for panel in panels )
y_min, y_max = min( panel.ya for panel in panels ), max( panel.ya for panel in panels )
x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)
y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)
size = 10
plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
plt.grid(True)
plt.xlabel('x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.plot(x, y, color='k', linestyle='-', linewidth=2)
plt.plot(np.append([panel.xa for panel in panels], panels[0].xa),
np.append([panel.ya for panel in panels], panels[0].ya),
linestyle='-', linewidth=1, marker='o', markersize=6, color='#CD2305');
Freestream conditions¶
The airfoil is immersed in a free-stream ($U_\infty$,$\alpha$) where $U_\infty$ and $\alpha$ are the velocity magnitude and angle of attack, respectively. Like before, we create a class for the free stream, even though we will only have one object that uses this class. It makes it easier to pass the free stream to other functions later on.
class Freestream:
"""Freestream conditions."""
def __init__(self, u_inf=1.0, alpha=0.0):
"""Sets the freestream conditions.
Arguments
---------
u_inf -- Farfield speed (default 1.0).
alpha -- Angle of attack in degrees (default 0.0).
"""
self.u_inf = u_inf
self.alpha = alpha*pi/180 # degrees --> radians
# defines and creates the object freestream
u_inf = 1.0
alpha = 0.0
freestream = Freestream(u_inf, alpha)
Flow tangency condition¶
A constant vortex strength $\gamma$ will be added to each panel (all panels have the same, constant vortex-sheet strength). Thus, using the principle of superposition, the velocity potential becomes:
$$ \begin{align*} \phi\left(x_{c_i},y_{c_i}\right) &= V_\infty x_{c_i} \cos \alpha + V_\infty y_{c_i} \sin \alpha \\ &+ \sum_{j=1}^N \frac{\sigma_j}{2\pi} \int_j \ln \left(\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \right) {\rm d}s_j \\ &- \sum_{j=1}^N \frac{\gamma}{2\pi} \int_j \tan^{-1} \left(\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\right) {\rm d}s_j \end{align*} $$We already worked the first integral in the previous lesson:
$$\frac{\partial}{\partial n_i} \ln \left(\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \right) = \frac{\left(x_{c_i}-x_j\right)\frac{\partial x_{c_i}}{\partial n_i} + \left(y_{c_i}-y_j\right)\frac{\partial y_{c_i}}{\partial n_i}}{\left(x_{c_i}-x_j\right)^2 + \left(x_{c_i}-x_j\right)^2}$$where $\frac{\partial x_{c_i}}{\partial n_i} = \cos \beta_i$ and $\frac{\partial y_{c_i}}{\partial n_i} = \sin \beta_i$, and:
$$x_j(s_j) = x_{b_j} - s_j \sin \beta_j$$$$y_j(s_j) = y_{b_j} + s_j \cos \beta_j$$We now need to derive the last integral of the boundary equation:
$$\frac{\partial}{\partial n_i} \tan^{-1} \left(\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\right)= \frac{\left(x_{c_i}-x_j\right)\frac{\partial y_{c_i}}{\partial n_i} - \left(y_{c_i}-y_j\right)\frac{\partial x_{c_i}}{\partial n_i}}{\left(x_{c_i}-x_j\right)^2 + \left(y_{c_i}-y_j\right)^2}$$where $\frac{\partial x_{c_i}}{\partial n_i} = \cos \beta_i$ and $\frac{\partial y_{c_i}}{\partial n_i} = \sin \beta_i$.
Kutta-condition¶
To enforce the Kutta-condition, we state that the pressure coefficient on the fisrt panel must be equal to that on the last panel:
$$C_{p_1} = C_{p_{N}}$$Using the definition of the pressure coefficient $C_p = 1-\left(\frac{V}{U_\infty}\right)^2$, the Kutta-condition implies that the magnitude of the velocity at the first panel center must equal the magnitude of the last panel center:
$$V_1^2 = V_N^2$$Since the flow tangency condition requires that $V_{n_1} = V_{n_N} = 0$, we end up with the following Kutta-condition:
$$V_{t_1} = - V_{t_N}$$(the minus sign comes from the reference axis we chose for the normal and tangential vectors).
Therefore, we need to evaluate the tangential velocity at the first and last panels. Let's derive it for every panel, since it will be useful to compute the pressure coefficient.
$$V_{t_i} = \frac{\partial}{\partial t_i} \left(\phi\left(x_{c_i},y_{c_i}\right)\right)$$i.e.,
$$ \begin{align*} V_{t_i} &= V_\infty \sin \left(\alpha-\beta_i\right) \\ &+ \sum_{j=1,j\neq i}^N \frac{\sigma_j}{2\pi} \int_j \frac{\partial}{\partial t_i} \ln \left(\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \right) {\rm d}s_j \\ &- \sum_{j=1,j\neq i}^N \frac{\gamma}{2\pi} \int_j \frac{\partial}{\partial t_i} \tan^{-1} \left(\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\right) {\rm d}s_j \end{align*} $$which gives
$$ \begin{align*} V_{t_i} &= V_\infty \sin \left(\alpha-\beta_i\right) \\ &+ \sum_{j=1,j\neq i}^N \frac{\sigma_j}{2\pi} \int_j \frac{\left(x_{c_i}-x_j\right)\frac{\partial x_{c_i}}{\partial t_i} + \left(y_{c_i}-y_j\right)\frac{\partial y_{c_i}}{\partial t_i}}{\left(x_{c_i}-x_j\right)^2 + \left(x_{c_i}-x_j\right)^2} {\rm d}s_j \\ &- \sum_{j=1,j\neq i}^N \frac{\gamma}{2\pi} \int_j \frac{\left(x_{c_i}-x_j\right)\frac{\partial y_{c_i}}{\partial t_i} - \left(y_{c_i}-y_j\right)\frac{\partial x_{c_i}}{\partial t_i}}{\left(x_{c_i}-x_j\right)^2 + \left(x_{c_i}-x_j\right)^2} {\rm d}s_j \end{align*} $$where $\frac{\partial x_{c_i}}{\partial t_i} = -\sin \beta_i$ and $\frac{\partial y_{c_i}}{\partial t_i} = \cos \beta_i$
Building the linear system¶
Here, we build and solve the linear system of equations of the form
$$[A][\sigma,\gamma] = [b]$$where the $N+1 \times N+1$ matrix $[A]$ contains three blocks: an $N \times N$ source matrix (the same one of Lesson 10), an $N \times 1$ vortex array to store the weight of the variable $\gamma$ at each panel, and a $1 \times N+1$ Kutta array that repesents our Kutta-condition.
We are going to re-use the function integral() from Lesson 10 to compute the different integrals with the NumPy function integrate.quad():
def integral(x, y, panel, dxdz, dydz):
"""Evaluates the contribution of a panel at one point.
Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.
dxdz -- derivative of x in the z-direction.
dydz -- derivative of y in the z-direction.
Returns
-------
Integral over the panel of the influence at one point.
"""
def func(s):
return ( ((x - (panel.xa-sin(panel.beta)*s))*dxdz
+ (y - (panel.ya+cos(panel.beta)*s))*dydz)
/ ((x - (panel.xa-sin(panel.beta)*s))**2
+ (y - (panel.ya+cos(panel.beta)*s))**2) )
return integrate.quad(lambda s:func(s), 0.0, panel.length)[0]
We first build the source matrix:
def source_matrix(panels):
"""Builds the source matrix.
Arguments
---------
panels -- array of panels.
Returns
-------
A -- NxN matrix (N is the number of panels).
"""
N = len(panels)
A = np.empty((N, N), dtype=float)
np.fill_diagonal(A, 0.5)
for i, p_i in enumerate(panels):
for j, p_j in enumerate(panels):
if i != j:
A[i,j] = 0.5/pi*integral(p_i.xc, p_i.yc, p_j, cos(p_i.beta), sin(p_i.beta))
return A
then the vortex array:
def vortex_array(panels):
"""Builds the vortex array.
Arguments
---------
panels - array of panels.
Returns
-------
a -- 1D array (Nx1, N is the number of panels).
"""
a = np.zeros(len(panels), dtype=float)
for i, p_i in enumerate(panels):
for j, p_j in enumerate(panels):
if i != j:
a[i] -= 0.5/pi*integral(p_i.xc, p_i.yc, p_j, +sin(p_i.beta), -cos(p_i.beta))
return a
and finally the Kutta array:
def kutta_array(panels):
"""Builds the Kutta-condition array.
Arguments
---------
panels -- array of panels.
Returns
-------
a -- 1D array (Nx1, N is the number of panels).
"""
N = len(panels)
a = np.empty(N+1, dtype=float)
a[N] = 0.0
a[0] = 0.5/pi*integral(panels[N-1].xc, panels[N-1].yc, panels[0],
-sin(panels[N-1].beta), +cos(panels[N-1].beta))
a[N-1] = 0.5/pi*integral(panels[0].xc, panels[0].yc, panels[N-1],
-sin(panels[0].beta), +cos(panels[0].beta))
for i, panel in enumerate(panels[1:N-1]):
a[i] = 0.5/pi*(integral(panels[0].xc, panels[0].yc, panel,
-sin(panels[0].beta), +cos(panels[0].beta))
+ integral(panels[N-1].xc, panels[N-1].yc, panel,
-sin(panels[N-1].beta), +cos(panels[N-1].beta)) )
a[N] -= 0.5/pi*(integral(panels[0].xc, panels[0].yc, panel,
+cos(panels[0].beta), +sin(panels[0].beta))
+ integral(panels[N-1].xc, panels[N-1].yc, panel,
+cos(panels[N-1].beta), +sin(panels[N-1].beta)) )
return a
Now that the three blocks have be defined, we can assemble the matrix $[A]$:
def build_matrix(panels):
"""Builds the matrix of the linear system.
Arguments
---------
panels -- array of panels.
Returns
-------
A -- (N+1)x(N+1) matrix (N is the number of panels).
"""
N = len(panels)
A = np.empty((N+1, N+1), dtype=float)
AS = source_matrix(panels)
av = vortex_array(panels)
ak = kutta_array(panels)
A[0:N,0:N], A[0:N,N], A[N,:] = AS[:,:], av[:], ak[:]
return A
On the right hand-side, we store the free-stream conditions that do not depend on the unknowns strengths:
def build_rhs(panels, freestream):
"""Builds the RHS of the linear system.
Arguments
---------
panels -- array of panels.
freestream -- farfield conditions.
Returns
-------
b -- 1D array ((N+1)x1, N is the number of panels).
"""
N = len(panels)
b = np.empty(N+1,dtype=float)
for i, panel in enumerate(panels):
b[i] = - freestream.u_inf * cos(freestream.alpha - panel.beta)
b[N] = -freestream.u_inf*( sin(freestream.alpha-panels[0].beta)
+sin(freestream.alpha-panels[N-1].beta) )
return b
All in all, the system has been defined with the functions source_matrix(), vortex_array(), kutta_array() and build_rhs(), which we have to call to build the complete system:
A = build_matrix(panels) # computes the singularity matrix
b = build_rhs(panels, freestream) # computes the RHS array
Solving the linear system¶
The linear system is then solved using the NumPy function linalg.solve()and we store the results in the attribute sigma of each object. We also create a variable gamma to store the value of the constant vortex strength.
# solves the linear system
variables = np.linalg.solve(A, b)
for i, panel in enumerate(panels):
panel.sigma = variables[i]
gamma = variables[-1]
Surface pressure coefficient¶
The pressure coefficient at the $i$-th panel center is:
$$C_{p_i} = 1 - \left(\frac{V_{t_i}}{U_\infty}\right)^2$$So, we have to compute the tangential velocity at each panel center using the function get_tangential_velocity():
def get_tangential_velocity(panels, freestream, gamma):
"""Computes the tangential velocity on the surface.
Arguments
---------
panels -- array of panels.
freestream -- farfield conditions.
gamma -- circulation density.
"""
N = len(panels)
A = np.empty((N, N+1), dtype=float)
np.fill_diagonal(A, 0.0)
A[:,N] = 0.0
for i, p_i in enumerate(panels):
for j, p_j in enumerate(panels):
if i != j:
A[i,j] = 0.5/pi*integral(p_i.xc, p_i.yc, p_j, -sin(p_i.beta), +cos(p_i.beta))
A[i,N] -= 0.5/pi*integral(p_i.xc, p_i.yc, p_j, +cos(p_i.beta), +sin(p_i.beta))
b = freestream.u_inf * np.sin([freestream.alpha - panel.beta for panel in panels])
var = np.append([panel.sigma for panel in panels], gamma)
vt = np.dot(A, var) + b
for i, panel in enumerate(panels):
panel.vt = vt[i]
# computes the tangential velocity at each panel center.
get_tangential_velocity(panels, freestream, gamma)
def get_pressure_coefficient(panels, freestream):
"""Computes the surface pressure coefficients.
Arguments
---------
panels -- array of panels.
freestream -- farfield conditions.
"""
for panel in panels:
panel.cp = 1.0 - (panel.vt/freestream.u_inf)**2
# computes surface pressure coefficients
get_pressure_coefficient(panels, freestream)
# plots the surface pressure coefficients
val_x, val_y = 0.1, 0.2
x_min, x_max = min(panel.xa for panel in panels), max(panel.xa for panel in panels)
cp_min, cp_max = min(panel.cp for panel in panels), max(panel.cp for panel in panels)
x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)
y_start, y_end = cp_min-val_y*(cp_max-cp_min), cp_max+val_y*(cp_max-cp_min)
plt.figure(figsize=(10, 6))
plt.grid(True)
plt.xlabel('x', fontsize=16)
plt.ylabel('$C_p$', fontsize=16)
plt.plot([panel.xc for panel in panels if panel.side == 'extrados'],
[panel.cp for panel in panels if panel.side == 'extrados'],
color='r', linestyle='-', linewidth=2, marker='o', markersize=6)
plt.plot([panel.xc for panel in panels if panel.side == 'intrados'],
[panel.cp for panel in panels if panel.side == 'intrados'],
color='b', linestyle='-', linewidth=1, marker='o', markersize=6)
plt.legend(['extrados', 'intrados'], loc='best', prop={'size':14})
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.gca().invert_yaxis()
plt.title('Number of panels : %d' % N);
Accuracy check¶
For a closed body, the sum of all the source strengths must be zero. If not, it means the body would be adding or absorbing mass from the flow! Therfore, we should have
$$\sum_{i=1}^{N} \sigma_i l_i = 0$$where $l_i$ is the length of the $i^{\text{th}}$ panel.
With this, we can get a measure of the accuracy of the source panel method.
# calculates the accuracy
accuracy = sum(panel.sigma*panel.length for panel in panels)
print '--> sum of source/sink strengths:', accuracy
--> sum of source/sink strengths: 0.00197648461643
Lift coefficient¶
The lift is given by the Kutta-Joukowski theorem, $L = \rho \Gamma U_\infty$, where $\rho$ is the fluid density. The total circulation $\Gamma$ is given by:
$$\Gamma = \sum_{i=1}^N \gamma l_i$$Finally, the lift coefficient is given by:
$$C_l = \frac{\sum_{i=1}^N \gamma l_i}{\frac{1}{2}U_\infty c}$$with $c$ the chord-length of the airoil
# calculates of the lift
cl = gamma*sum(panel.length for panel in panels)/(0.5*freestream.u_inf*(x_max-x_min))
print '--> Lift coefficient: Cl = %.3f' % cl
--> Lift coefficient: Cl = 0.008
Boundary Layer Calculation¶
To calculate the boundary layer in the laminar region, we use Thwaites' method. The laminar/tubulent transition is prediction with Michel's criterion. Finally, the turbulent boundary layer is solved using Head's method.
We start by finding the stagnation point on the surface near the leading edge with the function find_stagnation(). The function return the coordinates of the stagnation point (x_stagn, y_stagn) and the index of the station I just after the stagnation point.
def find_stagnation(panels):
"""Finds the stagnation point.
Arguments
---------
panels -- array of panels.
Returns
-------
I -- index of the pannel control-point just after the stagnation point.
x_stagn, y_stagn -- Cartesian coordinates of the stagnation point.
"""
I = 0
for i, panel in enumerate(panels):
if panel.vt/panels[0].vt < 0.0:
I = i
break
# station before stagnation point
vt1, x1, y1 = panels[I-1].vt, panels[I-1].xc, panels[I-1].yc
# station after stagnation point
vt2, x2, y2 = panels[I].vt, panels[I].xc, panels[I].yc
# interpolation to find stagnation point
x_stagn, y_stagn = x1-vt1*(x2-x1)/(vt2-vt1), y1-vt1*(y2-y1)/(vt2-vt1)
return I, x_stagn, y_stagn
# computes the stagnation point
I, x_stagn, y_stagn = find_stagnation(panels)
# plots the geometry and the stagnation point
val_x, val_y = 0.1, 0.2
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()
x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)
y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)
size = 10
plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
plt.grid(True)
plt.xlabel('x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.plot(x, y, color='k', linestyle='-', linewidth=1)
plt.scatter(x_stagn, y_stagn, s=40, color='r', marker='o', linewidth=0)
plt.title('Stagnation point', fontsize=16);
Description of the algorithm implemented¶
The boundary layer calculation has been embeded in a single class BoundaryLayer. The constructor of the class builds the running coordinate s along a boundary, stores the edge velocity from the source-vortex panel method described above and calculates the velocity gradient.
For each instance of the class BoundaryLayer, we start by calculating the momentum thickness theta, the shape factor H and the friction coefficient cf by integration along the boundary, marching from the stagnation point (using a 4th-order Runge-Kutta method).
The process starts with Thwaites' method (class' method solve_laminar()). After each station, Michel's criterion is applied to check is laminar/turbulent transition happened (class' method transition_prediction()). If we are still in the laminar region, Thwaites' method is applied to the next station. If we are now in the turbulent region, Head's method is used to continue the integration along the boundary (class' method solve_turbulent()).
After having solved for the momentum thickness, the shape factor and the friction coefficient, the displacement thickness is calculated and plotted (class' method displacement_thickness()).
Thwaites' Method for Laminar Boundary Layer¶
Thwaites' method relies on the momentum integral equation but makes no assumption on the form of the velocity profile.
The method supplements the momentum integral equation with algebraic relations among the unknowns $\theta$ (momentum thickness), $H$ (shape factor) and $c_f$ (skin friction coefficient).
The method introduces the parameter $l$:
$$l \equiv \frac{1}{2} Re_{\theta} c_f$$where $Re_{\theta} = \frac{u_e \theta}{\nu}$ is the Reynolds number based on the momentum thickness.
Since the skin friction coefficient $c_f$ is a strong function of the Reynolds number, this new paramter $l$ becomes less dependent of the Reynolds number.
The momentum integral equation is multiplied by $Re_{\theta}$ to get:
$$Re_{\theta} \frac{d\theta}{ds} + Re_{\theta} (2+H) \frac{\theta}{u_e} \frac{du_e}{ds} = l$$Thwiates defined a dimensionless pressure-gradient parameter:
$$\lambda \equiv \frac{\theta^2}{\nu} \frac{du_e}{ds}$$so that the modified momentum integral relation becomes:
$$\frac{u_e}{\nu}\frac{du_e^2}{ds} = 2 \left( l - (2+H) \lambda \right)$$Thwaites approximated the right hand-side with the linear formula:
$$2 \left( l - (2+H) \lambda \right) \approx 0.45 - 6\lambda$$that, when used, leads to the following equation:
$$\frac{d}{ds} \left( \theta^2 u_e^6 \right) = 0.45 \frac{u_e^5}{\nu}$$The above equation can be numerically solved to get the momentum thickness (in the laminar region), using a 4th-order Runge-Kutta method with the given initial condition:
$$ \theta(s=0) = \sqrt{\frac{0.075\nu}{\frac{du_e}{ds}(0)}} $$Then, the pressure-gradient parameter $\lambda$ can be calculated.
Finally, the shape factor $H$ and the skin friction coefficient $c_f = \frac{2l}{Re_{\theta}}$ are calculated using correlations from Cebecci and Bradshaw:
$$ H(\lambda) = 2.61 - 3.75 \lambda + 5.24 \lambda^2 \quad\quad \text{if} \quad 0 < \lambda < 0.1 $$$$ H(\lambda) = 2.088 + \frac{0.0731}{\lambda+0.14} \quad\quad \text{if} \quad -0.1 < \lambda < 0 $$and
$$ l(\lambda) = 0.22 + 1.57\lambda -1.8 \lambda^2 \quad\quad \text{if} \quad 0 < \lambda < 0.1 $$$$ l(\lambda) = 0.22 + 1.402\lambda +\frac{0.018\lambda}{\lambda+ 0.107} \quad\quad \text{if}\quad -0.1 < \lambda < 0 $$Michel's Transition Criterion¶
The laminar/turbulent transition point is determined using Michel's criterion which states that transition occurs in the boundary layer when:
$$ Re_{\theta} \geq 1.174 \left( Re_s^{0.46} + 22,400 Re_s^{-0.54} \right) $$where $Re_s = \frac{u_e s}{\nu}$ is the Reynolds number based on the position along the airfoil ($s$ is the running coordinate).
Head's Entrainment Method for Turbulent Boundary Layers¶
To solve the turbulent boundary-layer, we use the Von Karman momentum integral relation:
$$ \frac{d\theta}{ds} + (2+H) \frac{\theta}{u_e} \frac{du_e}{ds} = \frac{1}{2} c_f $$where the skin fricition coefficient is approximated using Ludwieg-Tillmann relation:
$$ c_f \approx 0.246 Re_{\theta}^{-0.268} 10^{-0.678H} $$The entrainment relation is used:
$$ E = \frac{1}{u_e} \frac{d}{ds} \left( u_e (\delta-\delta^*) \right) $$$E$ needs to be correlated with other integral parameters. The original method of Head (1958) defines a new shape factor:
$$ H_1 =\frac{\delta-\delta^*}{\theta} $$and writes the entrainment equation as a curve-fit correlation:
$$ \frac{d}{ds}\left( u_e \theta H_1 \right) = u_e F(H_1) $$where
$$ F(H_1) \approx 0.0306 \left( H_1 - 3.0 \right)^{-0.6169} $$Head showed that $H_1$ is correlated to the standard shape factor $H$:
$$ H_1 = 3.3 + 0.8234(H-1.1)^{-1.287} \quad\quad \text{for} \quad H \leq 1.6 $$$$ H_1 = 3.3 + 1.5501(H-0.6678)^{-3.064} \quad\quad \text{for} \quad H \geq 1.6 $$Using the momentum integral relation, the entrainment can be written as follow:
$$ \frac{dH_1}{ds} = (1+H)\frac{H_1}{u_e}\frac{du_e}{ds} + \frac{1}{\theta}F(H_1) - \frac{H_1}{\theta}\frac{c_f}{2} $$The above equation and the momentum integral relation are numericall y solved in the turbulent region using a 4th-order Runge-Kutta method.
class BoundaryLayer:
"""Boundary-layer correction."""
def __init__(self, x_stagn, y_stagn, panels, i_start, i_end):
"""Initializes the boundary layer,
computes the running coordinate, the edge-velocity and its gradient.
Arguments
---------
x_stagn, y_stagn -- Cartesian coordinates of the stagnation point.
panels -- array of panels.
i_start, i_end -- indices of the first and last panels on the boundary.
"""
self.N = abs(i_end-i_start) + 2
# stores panels information related to a side
# the parameter 'a' is used to reverse the sign of the tangential velocity
# due to sign convention in panel method
if i_start > i_end:
self.panels = np.flipud(panels[i_end:i_start+1])
a = -1.
else:
self.panels = panels[i_start:i_end+1]
a = +1.
# stores the Cartesian coordinates
self.x = np.append(x_stagn, [panel.xc for panel in self.panels])
self.y = np.append(y_stagn, [panel.yc for panel in self.panels])
self.beta = np.append(pi, [panel.beta for panel in self.panels])
# computes running coordinate
self.s = np.empty(self.N, dtype=float)
self.s[0] = 0.0
self.s[1] = sqrt( (x_stagn-self.panels[0].xc)**2 + (y_stagn-self.panels[0].yc)**2 )
for i in range(2,self.N):
self.s[i] = self.s[i-1] + self.panels[i-2].length/2 + self.panels[i-1].length/2
# stores the tangetial surface velocity and computes the gradient
self.ve = np.append(0.0, [a*panel.vt for panel in self.panels])
self.velocity_gradient()
def solve(self, nu=1.0E-05):
"""Solves the boundary layer.
Arguments
---------
nu -- kinematic viscosity of the fluid (default 1.0E-05).
"""
self.nu = nu
# initializes the boundary layer variables
self.theta = np.zeros(self.N, dtype=float)
self.H = np.zeros(self.N, dtype=float)
self.cf = np.zeros(self.N, dtype=float)
# solves laminar boundary
self.solve_laminar()
# solves turbulent boundary
if not self.is_laminar:
self.solve_turbulent()
def solve_laminar(self):
"""Solves the laminar boundary layer (Thwaites' method)."""
def shape_factor(lamBda):
"""Returns the shape factor (Cebeci and Bradshaw correlations).
Arguments
---------
lamBda -- Thwaites pressure-gradient parameter.
"""
if 0. <= lamBda < 0.1:
return 2.61 - 3.75*lamBda + 5.24*lamBda**2
elif -0.1 < lamBda < 0.:
return 2.088 + 0.0731/(lamBda+0.14)
def skin_friction(lamBda, Re_theta):
"""Returns the skin friction (Cebeci and Bradshaw correlations).
Arguments
---------
lamBda -- Thwaites pressure-gradient parameter.
Re_theta -- Reynolds number based on momentum thickness.
"""
if 0.0 <= lamBda < 0.1:
l = 0.22 + 1.57*lamBda - 0.18*lamBda**2
elif -0.1 < lamBda < 0.0:
l = 0.22 + 1.402*lamBda + 0.018*lamBda/(lamBda+0.107)
return 2*l/Re_theta
def f(s, y):
"""Returns right hand-side of first-order ODE.
Arguments
---------
s -- running coordinate along the boundary.
y -- variable to be integrated.
"""
ve = (s-self.s[i-1])/h*(self.ve[i]-self.ve[i-1]) + self.ve[i-1]
dveds = (s-self.s[i-1])/h*(self.dveds[i]-self.dveds[i-1]) + self.dveds[i-1]
rhs = 0.45*self.nu/ve - 6*y/ve*dveds
return rhs
# initial conditions
self.theta[0] = sqrt(0.075*self.nu/abs(self.dveds[0]))
self.ve[0] = 1.0E-06
y = np.array([self.theta[0]**2])
# marching along the surface with RK4 integration
self.is_laminar = True
i = 1
while self.is_laminar and i < self.N:
s = self.s[i-1]
h = self.s[i] - self.s[i-1]
k1 = h * f(s, y)
k2 = h * f(s + h/2, y + k1/2)
k3 = h * f(s + h/2, y + k2/2)
k4 = h * f(s + h, y + k3)
y += 1./6 * (k1 + 2*k2 + 2*k3 + k4)
#self.theta[i] = sqrt(y[0])
self.theta[i] = sqrt( 0.45*self.nu/self.ve[i]**6
* np.trapz(self.ve[:i]**5,self.s[:i]) )
# checks for transition laminar/turbulent
self.transition_prediction(i)
if self.is_laminar:
# calculates the pressure-gradient parameter
lamBda = self.theta[i]**2 / self.nu * self.dveds[i]
# computes the shape factor
self.H[i] = shape_factor(lamBda)
# calculates the Reynolds number based on momentum thickness
Re_theta = self.ve[i] * self.theta[i]/self.nu
# computes the friction coefficient
self.cf[i] = skin_friction(lamBda, Re_theta)
i += 1
def solve_turbulent(self):
"""Solves the turbulent boundary layer (Head's method)."""
def skin_friction(theta, H, ve):
"""Returns the skin friction (Ludwieg-Tillmann relation).
Arguments
---------
theta -- momentum thickness.
H -- shape factor.
ve -- edge velocity.
"""
return 0.246 * (ve*theta/self.nu)**(-0.268) * 10.**(-0.678*H)
def new_shape_factor(H):
"""Returns Head's shape factor (Head's curve-fit formula).
Arguments
---------
H -- shape factor.
"""
if H < 1.1:
return 16.0
elif H <= 1.6:
return 3.3 + 0.8234*(H-1.1)**(-1.287)
else:
return 3.3 + 1.5501*(H-0.6778)**(-3.064)
def shape_factor(H1):
"""Returns the shape factor (Head's curve-fit formula).
Arguments
---------
H1 -- Head's shape factor.
"""
if H1 <= 3.32:
return 3.0
elif H1 < 5.3:
return 0.6778 + ((H1-3.3)/1.5501)**(-(1./3.064))
else:
return 1.1 + ((H1-3.3)/0.8234)**(-(1./1.287))
def F(H1):
"""Returns RHS of entrainment equation (Head's curve-fit formula).
Arguments
---------
H1 -- Head's shape factor.
"""
return 0.0306 * (H1-3.0)**(-0.6169)
def f(s, y):
"""Returns right hand-side of first-order ODE.
Arguments
---------
s -- running coordinate along the boundary.
y -- variable to be integrated.
"""
theta, H1 = y[0], y[1]
ve = (s-self.s[i-1])/h*(self.ve[i]-self.ve[i-1]) + self.ve[i-1]
dveds = (s-self.s[i-1])/h*(self.dveds[i]-self.dveds[i-1]) + self.dveds[i-1]
H = shape_factor(H1)
cf = skin_friction(theta, H, ve)
rhs1 = 0.5*cf - (2.0 + H)*theta/ve*dveds
rhs2 = - H1/theta*rhs1 - H1/ve*dveds + F(H1)/theta
return np.array([rhs1, rhs2])
# initial conditions
H1 = new_shape_factor(self.H[self.i_trans-1])
y = np.array([self.theta[self.i_trans-1], H1])
# marching along the surface from the first turbulent station using RK4
for i in range(self.i_trans, self.N):
s = self.s[i-1]
h = self.s[i] - self.s[i-1]
k1 = h * f(s, y)
k2 = h * f(s + h/2, y + k1/2)
k3 = h* f(s + h/2, y + k2/2)
k4 = h * f(s + h, y + k3)
y += 1./6 * (k1 + 2*k2 + 2*k3 + k4)
self.theta[i], H1 = y[0], y[1]
self.H[i] = shape_factor(H1)
self.cf[i] = skin_friction(self.theta[i], self.H[i], self.ve[i])
def velocity_gradient(self):
"""Computes the velocity gradient w.r.t. the running coordinate."""
N = len(self.s)
self.dveds = np.gradient(self.ve,
np.append(self.s[1] - self.s[0],
self.s[1:N] - self.s[0:N-1]))
def transition_prediction(self, i):
"""Checks laminar/turbulent transition (Michel's criterion).
Arguments
---------
i -- index of the boundary point.
"""
Re_theta = self.ve[i]*self.theta[i]/self.nu
Re_s = self.ve[i]*self.s[i]/self.nu
michel_criterion = 1.174*(Re_s**(0.46)+22400*Re_s**(-0.54))
if Re_theta > michel_criterion:
self.is_laminar = False
self.i_trans = i # index of the first turbulent boundary point
def displacement_thickness(self, factor=1.0):
"""Computes the displacement thickness and its Cartesian coordinates.
Arguments
---------
factor -- amplification of the boundary thickness (default 1.0).
"""
self.delta_star = self.H * self.theta
self.x_delta = self.x + factor*self.delta_star*np.cos(self.beta)
self.y_delta = self.y + factor*self.delta_star*np.sin(self.beta)
Now that the stagnation has been found, the airfoil surface can be divided into two boundaries: upper and lower and we create two objects of the class BoundaryLayer:
# creates the different boundaries to solve
boundaries = {}
boundaries['upper'] = BoundaryLayer(x_stagn, y_stagn, panels, I-1, 0)
boundaries['lower'] = BoundaryLayer(x_stagn, y_stagn, panels, I, N-1)
For each surface, the boundary layer is calculated using the algorithm described previously. The Reynolds number has been set to $Re = 10^6$.
The displacement thickness $\delta^*$ is then calculated using the definition of the shape factor:
$$ H \equiv \frac{\delta^*}{\theta} $$Note that we use a factor $20$ to make the boundary layer visible on the figure below.
# solves the boundary layer and computes the displacement thickness on each side
Re = 1.0E+06
for boundary in boundaries.itervalues():
boundary.solve(nu=1./Re)
boundary.displacement_thickness(factor=20)
-c:89: RuntimeWarning: divide by zero encountered in double_scalars
Finally, the displacement thickness is displayed on the airfoil geometry. The solid blue line represents the boundary layer thickness, and the red dot is the laminar/turbulent transition point.
# plots the displacement thickness of the different boundaries
val_x, val_y = 0.1, 1.0
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()
x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)
y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)
size = 10
plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
plt.grid(True)
plt.xlabel('x', fontsize=16)
plt.ylabel('y', fontsize=16)
plt.xlim(x_start, x_end)
plt.ylim(y_start, y_end)
plt.plot(x, y, color='k', linestyle='-', linewidth=1)
for boundary in boundaries.itervalues():
plt.plot(boundary.x_delta, boundary.y_delta, color='b', linestyle='-', linewidth=1);
if not boundary.is_laminar:
plt.scatter(boundary.x_delta[boundary.i_trans],
boundary.y_delta[boundary.i_trans],
s=40, color='r', marker='o', zorder=10);
Example of the flat plate¶
To test the numerical method implemented above, we use the example of the flat plate to compare our data with experimental ones.
A 5-meter flat plate is immersed in a uniform flow (with speed $33$ m/s and viscosity of $1.51 \times 10^{-5}$ m2/s).
class Plate:
"""Contains plate information."""
def __init__(self, x_min, x_max, ve=1.0, nu=1.51E-05, N=100):
"""Discretizes the plate into panels and stores edge-velocity.
Arguments
---------
x_min, x_max -- boundaries of the plate.
ve -- edge-velocity (default 1.0).
nu -- kinematic viscosity (default 1.51E-05).
N - number of panels (default 100).
"""
self.x_min = x_min
self.x_max = x_max
self.N = N
self.x = np.linspace(self.x_min, self.x_max, self.N+1)
self.panels = np.empty(self.N, dtype=object)
for i, panel in enumerate(self.panels):
self.panels[i] = Panel(self.x[i], 0.0, self.x[i+1], 0.0)
self.panels[i].vt = ve
self.nu = nu
self.boundary = None
The flat plate is discretized into N=200 panels:
# creates the flat plate
plate = Plate(0.0, 5.0, nu=1.51E-05, ve=33.0, N=200)
And we create an object plate of the class BoundaryLayer to solve the boundary layer equations.
# defines the boundary layer
plate.boundary = BoundaryLayer(plate.x[0], 0.0, plate.panels, 0, plate.N-1)
# solves the boundary layer
plate.boundary.solve(nu=plate.nu)
The Blasius solution is also computed for comparison:
# Blasius solution
bound = plate.boundary
momentum_thickness = np.empty(bound.N, dtype=float)
momentum_thickness[0] = 0.0
if boundary.is_laminar:
momentum_thickness[1:bound.N] = 0.664*bound.s[1:bound.N] \
/ np.sqrt(bound.ve[1:bound.N]*bound.s[1:bound.N]
/ bound.nu)
else:
i_trans = bound.i_trans
momentum_thickness[1:i_trans] = 0.664 * bound.s[1:i_trans] \
/ np.sqrt(bound.ve[1:i_trans]*bound.s[1:i_trans]/bound.nu)
momentum_thickness[i_trans:bound.N] = 0.0155*bound.s[i_trans:bound.N] \
/ (bound.ve[i_trans:bound.N]
*bound.s[i_trans:bound.N]
/bound.nu)**(1./7)
Numerical results are compared with experimental data from Wieghardt and Tillmann (1951):
# reads experimental data from Wieghardt and Tillmann (1951)
with open ('../resources/WIECF.dat') as file_name:
data = np.loadtxt(file_name, dtype=float, usecols=(0, 3), skiprows=5)
x_exp = 0.3048*data[:,0] # conversion feet to meters
theta_exp = 0.3048*data[:,1] # conversion feet to meters
# plots the momentum thickness along the flat plate
plt.figure(figsize=(10, 6))
plt.grid(True)
plt.xlabel('x', fontsize=18)
plt.ylabel(r'$Re_\theta$', fontsize=18)
# plots numerical solution
plt.plot(bound.x, bound.theta*bound.ve/bound.nu,
color='b', linestyle='-', linewidth=1)
# plots Blasius solution
plt.plot(bound.s, momentum_thickness*bound.ve/bound.nu,
color='r', linestyle='-', linewidth=1)
# plots experimental data
plt.scatter(x_exp, theta_exp*bound.ve[1]/bound.nu, color='g', s=40, marker='o', linewidth=0)
plt.legend(['Present work', 'Blasius', 'Wieghardt & Tillmann (1951)'],
loc='best', prop={'size':16})
if not plate.boundary.is_laminar:
plt.axvline(bound.s[plate.boundary.i_trans], ymin=0., ymax=1.,
color='k', linestyle='-', linewidth=2)
The cell below can be ignored. It is used to create the style of the notebook.
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom.css', 'r').read()
return HTML(styles)
css_styling()