Method of characteristics¶
The method of characteristics is a numerical approach to modeling two-dimensional supersonic flow. (It can also be extended to axisymmetric and three-dimensional flows.)
import numpy as np
from scipy.optimize import root_scalar
import matplotlib.pyplot as plt
# these lines are only for helping improve the display
import matplotlib_inline.backend_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('pdf', 'png')
plt.rcParams['figure.dpi']= 400
plt.rcParams['savefig.dpi'] = 400
plt.rcParams['mathtext.fontset'] = 'cm'
We start with the two-dimensional, steady, velocity potential equation, which applies to irrotational and adiabatic flows:
where \(\Phi = f(x,y)\) is the velocity potential function, subscripts of \(x\) and \(y\) represent partial derivatives, and \(a\) is the local speed of sound. The velocity potential is related to velocity via
where the velocity vector is \( \mathbf{V} (x,y) = u \mathbf{\hat{i}} + v \mathbf{\hat{j}} \).
Since the partial derivatives of the velocity potential are also functions of position—in other words, \(\Phi_x = f(x,y)\) and \(\Phi_y = g(x,y)\)—we can take exact differentials of them:
Expressing these using velocity components, we can simplify to
If we then rewrite Equation (87) using the velocity components:
Now, with Equations (88), (89), and (90), we have a complete linear system of equations for solving \( \mathbf{\Phi} = \left[ \Phi_{xx}, \Phi_{xy}, \Phi_{yy} \right] \).
These unknowns are solveable using Cramer’s rule! For example, we can solve for \(\Phi_{xy}\):
At a particular location in space, \(\Phi_{xy}\) has a specific value, and general choices of \(dx\) and \(dy\) point away from this location. Along the line formed by \(dx\) and \(dy\), the velocity changes with particular \(du\) and \(dv\). Regardless of the choice of this line, the associated values of change in velocity will ensure that \(\Phi_{xy}\) remains constant.
Fig. 11 Geometry of a streamline.¶
Figure 11 shows the geometry of a streamline, with \(dx\) and \(dy\) directed away from point A, as well as the initial velocity \(V\) and components \(du\) and \(dv\).
Characteristic lines¶
However, a particular line in space, given by particular values of \(dx\) and \(dy\), leads to \(D = 0\). Since we know that \(\Phi_{xy}\) cannot physically be infinite at this location, then \(N = 0\) to keep \(\Phi_{xy}\) finite:
Or, the velocity derivatives are indeterminate in the directions created by \(dx\) and \(dy\). These lines in space are the characteristic lines, where the derivatives of velocity are indeterminate. We now have equations we can use to create these lines.
Setting \(D = 0\) leads to
where \( \left( \frac{dy}{dx} \right)_{\text{char}}\) gives the slope of the characteristic line. We can solve this quadratic equation to get
Equation (91) describes the slope of the characteristic line in space; if we focus on supersonic flows (where \( \frac{u^2 + v^2}{a^2} = M^2 > 1 \)), there are two real solutions to this equation.
Since \( u = V \cos \theta \) and \( v = V \sin \theta \), we can simplify the slope to
In other words, there are two characteristic lines passing through any point, which are also Mach lines. The left-running \( \theta + \mu \) line is the \(C_{+}\) characteristic, while the right-running \( \theta - \mu \) line is the \(C_{-}\) characteristic.
Compatibility equations¶
We now can get compatibility equations for finding how the flow velocity magnitude and direction relate along the characteristic lines.
Along a characteristic line, where \(D = 0\), we can set \(N = 0\), to get
By recognizing that \(\frac{dy}{dx} = \left(\frac{dy}{dx}\right)_{\text{char}}\), and incorporating Equation (92), we can get
Equation (93) describes how the flow changes along characteristic lines, relating changes in direction (\(\theta\)) to changes in magnitude (\(V\)). Integrating this equation actually gives the Prandtl–Meyer function \(\nu(M)\), and we can obtain algebraic compatibility equations:
where \(K_{-}\) and \(K_{+}\) are constants along the \(C_{-}\) and \(C_{+}\) characteristic lines, respectively, and the Prandtl–Meyer function is
Unit processes¶
We now have the tools to apply the method of characteristics to calculate how supersonic flow varies in a particular domain. The ways we apply these equations in specific situations are called unit processes.
Internal flow¶
Fig. 12 Internal flow unit process, with two characteristic lines intersecting to form a third point.¶
Figure 12 shows two characteristic lines, passing through points 1 and 2 where we know the flow conditions, intersecting at point 3.
Since we know the conditions at point 1, we know \(\nu_1\) and \(\theta_1\), and so we know the constant that holds along the entire \(C_{-}\) line:
and similarly since we know the conditions at point 2, we know the constant that holds along the \(C_{+}\) line:
The constants associated with both characteristic lines apply at point 3, and so from Equations (94) and (95) we have
which we can solve for the unknown values at point 3:
Once we have \(\nu_3\), we can solve for \(M_3\) using the Prandtl–Meyer function, and from that determine other properties such as pressure, temperature, and density via isentropic flow relations.
To actually determine the location of point 3, we can approximate that the characteristic lines are straight lines between the two points, and use the average slope based on the angles at the points. In other words, the \(C_{-}\) characteristic line connecting points 1 and 3 has the average slope
and the \(C_{+}\) characteristic line connecting points 2 and 3 has the average slope
Using the known coordinates of points 1 and 2, we can use the intersection of the two characteristic lines to find the location of point 3:
which we can solve to get
Example¶
Consider two points in a supersonic flow of air, both with Mach number of 2 and right above each other (i.e., \(\Delta x = 0\) between the points). Point 2 is along a line of symmetry and so the flow direction is horizontal. Point 1 is 0.25 above, and \(\theta_1 = 4°\). Find the properties at and location of point 3.
def get_mach_angle(mach):
'''Returns Mach angle'''
return np.arcsin(1.0 / mach)
def get_prandtl_meyer(mach, gamma=1.4):
'''Evaluate Prandtl-Meyer function at given Mach number.
Defined as the angle from the flow direction where Mach = 1 through which the
flow turned isentropically reaches the specified Mach number.
'''
return (
np.sqrt((gamma + 1) / (gamma - 1)) *
np.arctan(np.sqrt((gamma - 1)*(mach**2 - 1)/(gamma + 1))) -
np.arctan(np.sqrt(mach**2 - 1))
)
def solve_prandtl_meyer(mach, nu, gamma=1.4):
'''Solve for unknown Mach number, given Prandtl-Meyer function (in radians).'''
return (nu - get_prandtl_meyer(mach, gamma))
def get_reference_area(mach, gamma=1.4):
'''Calculate reference area ratio'''
return ((1.0/mach) * ((1 + 0.5*(gamma-1)*mach**2) /
((gamma + 1)/2))**((gamma+1) / (2*(gamma-1)))
)
def solve_mach_area(mach, area_ratio, gamma=1.4):
'''Used to find Mach number for given reference area ratio and gamma'''
return (area_ratio - get_reference_area(mach, gamma))
gamma = 1.4
M_1 = 2.0
theta_1 = 4 * np.pi/180
x1 = 0
y1 = 0.25
M_2 = 2.0
theta_2 = 0
x2 = 0
y2 = 0
nu_1 = get_prandtl_meyer(M_1, gamma)
nu_2 = get_prandtl_meyer(M_2, gamma)
K_minus = theta_1 + nu_1
K_plus = theta_2 - nu_2
theta_3 = 0.5 * (K_minus + K_plus)
nu_3 = 0.5 * (K_minus - K_plus)
root = root_scalar(solve_prandtl_meyer, x0=2.0, x1=2.5, args=(nu_3, gamma))
M_3 = root.root
print(f'M_3 = {M_3: .3f}')
print(f'θ_3 = {theta_3 * 180/np.pi: .3f}°')
dydx_minus = np.tan(
0.5*(theta_1+theta_3) - 0.5*(get_mach_angle(M_1) + get_mach_angle(M_3))
)
dydx_plus = np.tan(
0.5*(theta_2+theta_3) + 0.5*(get_mach_angle(M_2) + get_mach_angle(M_3))
)
x3 = (y2 - y1 - x2*dydx_plus + x1*dydx_minus) / (dydx_minus - dydx_plus)
y3 = y1 + (x3 - x1) * dydx_minus
plt.subplots(figsize=(6, 2))
plt.plot([x1, x3], [y1, y3], '-o')
plt.plot([x2, x3], [y2, y3], '-or')
plt.xlim([-0.025, 0.25])
plt.ylim([-0.04, 0.29])
plt.grid(True)
plt.tight_layout()
plt.show()
M_3 = 2.073
θ_3 = 2.000°
We can see the characteristic lines and grid points where we have the flow properties.
In general, the characteristic lines will be curved, but we can approximate the lines as straight between grid points.
Wall point¶
Fig. 13 Wall point unit process, with one characteristic line intersecting with a point on a solid wall.¶
Figure 13 shows a characteristic line, passing through point 4 and intersecting with a solid wall at point 5.
Since the flow properties are known at point 4, we know \(K_{-}\) for the \(C_{-}\) characteristic line, which is the same at point 5:
Since we know the shape of the wall, and the flow must be tangent to the wall, we know \(\theta_5\), and so we can solve for \(\nu_5\):
Example¶
Below a diverging-channel wall inclined at 10°, the Mach number of a horizontal flow of air is 2.0. This grid point (1) is located 0.2 below the wall, as Figure 14 shows. Find the conditions at the wall, at grid point 2.
Fig. 14 Grid point and characteristic line intersecting a wall.¶
We know \(M_1 = 2\) and \(\theta_1 = 0°\), and \(\theta_2 = 10°\). Since the grid point is located below the wall, points 1 and 2 are connected by a \(C_{+}\) characteristic line, where
so we can solve for \(\nu_2\):
gamma = 1.4
M_1 = 2
theta_1 = 0
nu_1 = get_prandtl_meyer(M_1, gamma)
theta_2 = 10 * np.pi/180
nu_2 = nu_1 + theta_2 - theta_1
root = root_scalar(solve_prandtl_meyer, x0=2.0, x1=2.5, args=(nu_2, gamma))
M_2 = root.root
print(f'M_2 = {M_2: .3f}')
M_2 = 2.385
To find the location of wall grid point 2 relative to flow grid point 1, we can consider that there is a point 0 located directly above point 1. We know the slope of the wall, and therefore we can use geometry to determine the location of the new point:
These equations can be solved similarly to how we got Equations (99) and (100) before:
dydx_02 = np.tan(theta_2)
dydx_plus = np.tan(theta_1 + get_mach_angle(M_1))
y0 = 0.2
x0 = 0
y1 = 0
x1 = 0
x2 = (y0 - y1 - x0*dydx_02 + x1*dydx_plus) / (dydx_plus - dydx_02)
y2 = y0 + (x2 - x0) * dydx_02
plt.subplots(figsize=(6, 2))
# construct line for wall
xs = np.linspace(-0.1, 0.6, 50, endpoint=True)
ys = y0 + np.tan(theta_2) * (xs - x0)
plt.plot([x1, x2], [y1, y2], '-o')
plt.plot(xs, ys, 'k')
plt.grid(True)
plt.tight_layout()
plt.show()
Example: diverging channel¶
Consider a situation where uniform radial flow at Mach 2.0 enters a two-dimensional diverging channel with straight walls. Compute the variation of Mach number in this radial flow field, assuming isentropic, steady flow. The walls are inclined at a total angle of 12° and \(\gamma = 1.4\).
Solution: We can take advantage of the symmetry of this problem and only solve for half the flow field, with the line of symmetry being treated as a horizontal wall. Then, the physical wall moves upward at an angle of 6°.
We choose a number of initial data points where the Mach number is 2.0. Since the flow is radial, the flow direction varies smoothly from 0° to 6°. Once we have this initial data line, we can proceed using the unit processes given above, treating the line of symmetry as a horizontal wall.
To find the coordinates of these initial points, we consider that the initial data line is an arc with length \(R \Delta \theta\), where \(\Delta \theta\) is the 6° angle of the half-channel. The upper wall is given by a line following
and we can fix the initial point (on the wall) with coordinate \(y_0 = 1\). Then, \(x_0 = 1 / \tan \Delta \theta\), and the radius of this arc is given by
The remaining points on the initial data line can be found with \( x = R \cos \theta \) and \( y = R \sin \theta \), and the remaining grid points throughout the flow field can be found using the equations given above.
gamma = 1.4
wall_angle = 6 * np.pi / 180
# initial data line
M_initial = 2.0
num_initial = 4
num_columns = 5
num_points = (
num_initial * num_columns +
(num_initial - 1) * (num_columns - 1)
)
thetas = np.zeros(num_points)
machs = np.zeros(num_points)
nus = np.zeros(num_points)
ys = np.zeros(num_points)
xs = np.zeros(num_points)
kplus = np.zeros(num_initial - 1)
kminus = np.zeros(num_initial - 1)
# initial line
delta_theta = wall_angle / (num_initial - 1)
thetas[:num_initial] = [
wall_angle - delta_theta * idx for idx in range(num_initial)
]
machs[:num_initial] = M_initial
nus[:num_initial] = get_prandtl_meyer(machs[:num_initial], gamma)
# initial data line (arc)
ys[0] = 1.0
xs[0] = 1 / np.tan(wall_angle)
radius_initial = np.sqrt(xs[0]**2 + ys[0]**2)
xs[1:num_initial] = radius_initial * np.cos(thetas[1:num_initial])
ys[1:num_initial] = radius_initial * np.sin(thetas[1:num_initial])
# calculate characteristic
for idx in range(num_initial - 1):
kminus[idx] = thetas[idx] + get_prandtl_meyer(machs[idx], gamma)
kplus[idx] = thetas[idx+1] - get_prandtl_meyer(machs[idx+1], gamma)
# next column
for idx in range(num_initial - 1):
thetas[num_initial + idx] = 0.5 * (kminus[idx] + kplus[idx])
nus[num_initial + idx] = 0.5 * (kminus[idx] - kplus[idx])
root = root_scalar(
solve_prandtl_meyer, x0=2.0, x1=3.0,
args=(nus[num_initial + idx], gamma)
)
machs[num_initial + idx] = root.root
mu = get_mach_angle(machs[num_initial + idx])
dydx_minus = np.tan(0.5*(
thetas[idx] + thetas[num_initial + idx] -
get_mach_angle(machs[idx]) - mu
))
dydx_plus = np.tan(0.5*(
thetas[idx + 1] + thetas[num_initial + idx] +
get_mach_angle(machs[idx + 1]) + mu
))
xs[num_initial + idx] = (
ys[idx+1] - ys[idx] -
xs[idx+1] * dydx_plus +
xs[idx] * dydx_minus
) / (dydx_minus - dydx_plus)
ys[num_initial + idx] = (
ys[idx] +
(xs[num_initial + idx] - xs[idx])*dydx_minus
)
# remaining
for icol in range(1, num_columns):
idx_start = icol * (2*num_initial - 1)
for idx in range(num_initial):
# wall point
if idx == 0:
thetas[idx_start + idx] = wall_angle
nus[idx_start + idx] = (
nus[idx_start + idx - num_initial + 1] +
wall_angle - thetas[idx_start + idx - num_initial + 1]
)
root = root_scalar(
solve_prandtl_meyer, x0=2.0, x1=3.0,
args=(nus[idx_start + idx], gamma)
)
machs[idx_start + idx] = root.root
mu = get_mach_angle(machs[idx_start + idx])
dydx_plus = np.tan(0.5*(
thetas[idx_start + idx - num_initial + 1] +
thetas[idx_start + idx] +
get_mach_angle(machs[idx_start + idx - num_initial + 1]) + mu
))
x1 = xs[idx_start + idx - num_initial + 1]
y1 = ys[idx_start + idx - num_initial + 1]
x0 = x1
y0 = ys[0] + (x0 - xs[0]) * np.tan(wall_angle)
xs[idx_start + idx] = (
y0 - y1 - x0*np.tan(wall_angle) + x1*dydx_plus
) / (dydx_plus - np.tan(wall_angle))
ys[idx_start + idx] = (
y0 + (xs[idx_start + idx] - x0) * np.tan(wall_angle)
)
# point on line of symmetry
elif idx == num_initial - 1:
thetas[idx_start + idx] = 0
nus[idx_start + idx] = (
nus[idx_start + idx - num_initial] +
thetas[idx_start + idx - num_initial]
)
root = root_scalar(
solve_prandtl_meyer, x0=2.0, x1=3.0,
args=(nus[idx_start + idx], gamma)
)
machs[idx_start + idx] = root.root
mu = get_mach_angle(machs[idx_start + idx])
dydx_minus = np.tan(0.5*(
thetas[idx_start + idx - num_initial] +
thetas[idx_start + idx] -
get_mach_angle(machs[idx_start + idx - num_initial]) - mu
))
x1 = xs[idx_start + idx - num_initial]
y1 = ys[idx_start + idx - num_initial]
xs[idx_start + idx] = (
x1 - y1 / dydx_minus
)
ys[idx_start + idx] = 0
# interior point
else:
thetas[idx_start + idx] = 0.5 * (kminus[idx-1] + kplus[idx])
nus[idx_start + idx] = 0.5 * (kminus[idx-1] - kplus[idx])
root = root_scalar(
solve_prandtl_meyer, x0=2.0, x1=3.0,
args=(nus[idx_start + idx], gamma)
)
machs[idx_start + idx] = root.root
mu = get_mach_angle(machs[idx_start + idx])
dydx_minus = np.tan(0.5*(
thetas[idx_start + idx - num_initial] +
thetas[idx_start + idx] -
get_mach_angle(machs[idx_start + idx - num_initial]) - mu
))
dydx_plus = np.tan(0.5*(
thetas[idx_start + idx - num_initial + 1] +
thetas[idx_start + idx] +
get_mach_angle(machs[idx_start + idx - num_initial + 1]) + mu
))
xs[idx_start + idx] = (
ys[idx_start + idx - num_initial + 1] -
ys[idx_start + idx - num_initial] -
xs[idx_start + idx - num_initial + 1]*dydx_plus +
xs[idx_start + idx - num_initial]*dydx_minus
) / (dydx_minus - dydx_plus)
ys[idx_start + idx] = (
ys[idx_start + idx - num_initial] +
(xs[idx_start + idx] - xs[idx_start + idx - num_initial])*dydx_minus
)
# new characteristic constants
kminus[1:] = kminus[:-1]
kminus[0] = thetas[idx_start] + nus[idx_start]
kplus[:-1] = kplus[1:]
kplus[-1] = (
thetas[idx_start + num_initial - 1] -
nus[idx_start + num_initial - 1]
)
idx_start += num_initial
if icol < num_columns - 1:
for idx in range(num_initial - 1):
thetas[idx_start + idx] = 0.5 * (kminus[idx] + kplus[idx])
nus[idx_start + idx] = 0.5 * (kminus[idx] - kplus[idx])
root = root_scalar(
solve_prandtl_meyer, x0=2.0, x1=3.0,
args=(nus[idx_start + idx], gamma)
)
machs[idx_start + idx] = root.root
mu = get_mach_angle(machs[idx_start + idx])
dydx_minus = np.tan(0.5*(
thetas[idx_start + idx - num_initial] +
thetas[num_initial + idx] -
get_mach_angle(machs[idx_start + idx - num_initial]) - mu
))
dydx_plus = np.tan(0.5*(
thetas[idx_start + idx - num_initial + 1] +
thetas[idx_start + idx] +
get_mach_angle(machs[idx_start + idx - num_initial + 1]) + mu
))
xs[idx_start + idx] = (
ys[idx_start + idx - num_initial + 1] -
ys[idx_start + idx - num_initial] -
xs[idx_start + idx - num_initial + 1]*dydx_plus +
xs[idx_start + idx - num_initial]*dydx_minus
) / (dydx_minus - dydx_plus)
ys[idx_start + idx] = (
ys[idx_start + idx - num_initial] +
(xs[idx_start+idx] - xs[idx_start + idx - num_initial])*dydx_minus
)
We can examine the flow properties at each location, and visualize the characteristic grid network.
plt.plot(xs, ys, 'o')
plt.plot(
[xs[num_initial-1], xs[-1]],
[ys[num_initial-1], ys[-1]], '-.k'
)
plt.plot(
[xs[0], xs[-num_initial]], [ys[0], ys[-num_initial]], '-b',
)
for idx, (x,y) in enumerate(zip(xs, ys)):
plt.text(
x, y + 0.02, f'{idx: d}',
verticalalignment='bottom', horizontalalignment='center'
)
plt.ylim([-0.05, np.max(ys)+0.1])
plt.tight_layout()
plt.show()
print('id Mach angle (°)')
for idx, (mach, theta) in enumerate(zip(machs, thetas)):
print(f'{idx: d} {mach: .4f} {theta*180/np.pi: .3f}')
id Mach angle (°)
0 2.0000 6.000
1 2.0000 4.000
2 2.0000 2.000
3 2.0000 -0.000
4 2.0365 5.000
5 2.0365 3.000
6 2.0365 1.000
7 2.0733 6.000
8 2.0733 4.000
9 2.0733 2.000
10 2.0733 0.000
11 2.1106 5.000
12 2.1106 3.000
13 2.1106 1.000
14 2.1483 6.000
15 2.1483 4.000
16 2.1483 2.000
17 2.1483 0.000
18 2.1864 5.000
19 2.1864 3.000
20 2.1864 1.000
21 2.2251 6.000
22 2.2251 4.000
23 2.2251 2.000
24 2.2251 0.000
25 2.2642 5.000
26 2.2642 3.000
27 2.2642 1.000
28 2.3039 6.000
29 2.3039 4.000
30 2.3039 2.000
31 2.3039 0.000
We can compare the results of the two-dimensional method of characteristics solution against the exact solution determined using our one-dimensional methods, using the reference area ratio.
Based on the known Mach and gas property at the entrance, we can calculate
and then use area ratios at the exit:
where we are exploiting the radial flow condition.
area_ratio_in = get_reference_area(machs[num_initial-1], gamma)
# find Mach number from one-dimensional exact solution
area_ratio_exit = (xs[-1] / xs[num_initial-1]) * area_ratio_in
root = root_scalar(
solve_mach_area, x0=2.0, x1=2.5, args=(area_ratio_exit, gamma)
)
M_exact = root.root
print(f'M_exit (exact) = {M_exact: .5f}')
print(f'M_exit (exact) = {machs[-1]: .5f}')
print(f'Error: {100*np.abs(M_exact - machs[-1])/M_exact: .3e} %')
M_exit (exact) = 2.30388
M_exit (exact) = 2.30387
Error: 4.018e-04 %
The method of characteristics solution matches the exact solution quite closely, even though we used a fairly course mesh.