Code
# import relevant libraries
import plotly.graph_objects as go
import numpy as np
from scipy.optimize import root_scalar
from IPython.display import Markdown
import plotly.io as pio
pio.renderers.default = 'iframe'Supersonic flow with friction
Introduction
In this notebook we will study how friction changes the Mach number in a constant area duct. As we will see the coefficient of friction plays the same role as the change in area of a converging frictionless duct, meaning if the flow is originally subsonic the Mach number will increase and if it is originally supersonic the Mach number will decrease.
Problem Workthrough
It will be useful to consider the following parameters associated with the duct:
| Parameters | Values |
|---|---|
| Inlet Mach number, \({\color[rgb]{0.041893,0.355669,0.727621}M}_{1}\) | 0.5 |
| Axial location at \({\color[rgb]{0.041893,0.355669,0.727621}M}_{1}\), \(x_1\) | 2 |
| Surface skin friction coefficient, \(c_f\) | 0.0025 |
| Duct diameter, \(D\) | 0.1 |
| Heat capacity ratio, \(\gamma\) | 1.4 |
Other quantities used here are as per the lecture slides.
Code
# Define constants and known values
M_1 = 0.5
x_1 = 2
c_f = 0.0025
D = 0.1
g = 1.4If the pipe continues onto infinity at some point the flow will choke (reach M=1), this is useful because there is an isentropic flow relation between the Mach number and the distance to the choked position:
Derivation for this formula may be found in the Appendix section.
\[ \large \require{color}\frac{4c_{f}}{D}L = \frac{1 - {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2} + \frac{\gamma +1}{2\gamma} \ \textrm{ln}\left(\frac{\frac{\gamma +1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}{1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}\right) \]
This finds the maximum required distance for the flow to choke, then since we know the distance between \(\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}\) and \(\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2\) (where we want to know the Mach number) we can also find the distance it takes for the flow to choke fro \(\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2\). Then using the distance bwetween \(\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2\) and \(\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}=1\) we can solve the above equation for Mach number to find \(\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2\). However, the equation above does not have an analytical solution so a process such as curve fitting or root finding must be used, in this notebook root finding using scipy is used.
Code
# define a function that calculates how far it takes for the flow to choke given the initial mach number
def max_L(M):
L = (D/(4*c_f))*((1-M**2)/(1.4*M**2) + (g+1)/(2*g) * np.log(((M**2)*(g+1)/2)/(1+(M**2)*(g-1)/2)))
return L
# define the length to where the flow chokes and a position we want to find the Mach number
L_choke = max_L(M_1)
x_2 = x_1 + L_choke/2
x=np.linspace(0, x_1+L_choke+1, 100)Consider the diagram below that shows the change in the Mach number as a function of the axial location.
Code
fig = go.Figure()
fig.add_scatter(x=x, y=np.e**(-x) + 1, fill ='tonexty', fillcolor='rgba(133, 163, 201, 0.5)', line_color='blue')
fig.add_scatter(x=x, y=-np.e**(-x) - 1, fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='blue')
fig.add_vline(x=x_1, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}$")
fig.add_vline(x=x_1+L_choke, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}=1$")
fig.add_vline(x=x_2, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{2}$")
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()Now, we define the length from \({\color[rgb]{0.041893,0.355669,0.727621}M}_{2}\) to where the flow chokes.
Code
# define the length from M_2 to where the flow chokes
L_2 = L_choke + x_1 - x_2
# define function to find Mach roots
def mach_function(M, L):
return (D/(4*c_f))*((1-M**2)/(1.4*M**2) + (g+1)/(2*g) * np.log(((M**2)*(g+1)/2)/(1+(M**2)*(g-1)/2))) - L
# using scipy find the root where L=L_2
M_2_root = root_scalar(mach_function, args=(L_2), x0=M_1, x1=1)
M_2 = M_2_root.root
Markdown(fr'The Mach number at a distance ${round(x_2 - x_1, 3)}m$ from $M_1$ the Mach value is ${round(M_2, 3)}$')The Mach number at a distance \(5.345m\) from \(M_1\) the Mach value is \(0.589\)
Appendix: derivation
Detailed derivations for the formulas used above are given below. Students need not memorize the formulas, but rather understand how they are derived from first principles.
Code
fig = go.Figure()
fig.add_shape(type='circle', xref='x', yref='y', x0=0, y0=0, x1=1, y1=3)
fig.add_shape(type='circle', xref='x', yref='y', x0=3, y0=0, x1=4, y1=3, line_dash='dash')
fig.add_shape(type='line', x0=0.5, y0=3, x1=3.5, y1=3)
fig.add_shape(type='line', x0=0.5, y0=0, x1=3.5, y1=0)
fig.add_annotation(x=0.5, ax=-0.5, y=1.5, ay=1.5, xref='x', axref='x', yref='y', ayref='y', arrowhead=5,
text=r"$p+\rho v^2$")
fig.add_annotation(x=3.5, ax=5.5, y=1.5, ay=1.5, xref='x', axref='x', yref='y', ayref='y', arrowhead=5,
text=r"$p + \frac{dp}{dx}dx+\rho v^2 +\frac{d\rho v^2}{dx}dx$")
fig.add_annotation(x=1.5, ax=2.5, y=3.2, ay=3.2, xref='x', axref='x', yref='y', ayref='y', arrowhead=5, text=r"$\tau_w$")
fig.add_annotation(x=1.5, ax=2.5, y=-0.2, ay=-0.2, xref='x', axref='x', yref='y', ayref='y', arrowhead=5)
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()The shear force is constant across the surface area of the pipe spanned by a small section of fluid \(dx\) this area is \(2\pi R dx = \pi Ddx\) where \(D\) is the diameter of the pipe. The pressure and velocity terms only act on the inlet and exit areas which is given by \(\pi R^2 = \pi \frac{D^2}{4}\). Firstly, we must sum the forces:
\[ \large \require{color}F=-{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w\pi Dd x + \left( {\color[rgb]{0.315209,0.728565,0.037706}p} +{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 \right) \pi \frac{D^2}{4} -\left({\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}dx + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 + \frac{d {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx}dx\right)\pi \frac{D^2}{4} = 0 \]
\[ \large \Rightarrow -\pi \frac{D^2}{4}dx\left(\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx}\right) - {\color[rgb]{0.990308,0.800015,0.121270}\tau}_w\pi Ddx = 0 \]
\[ \large \Rightarrow \require{color}-{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w\pi D = \pi \frac{D^2}{4}\left(\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{d({\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2)}{dx}\right) \]
\[ \large \Rightarrow -\frac{4{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w}{D} = \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx} \]
If we assume the flow conditions
\[ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} = const, \; \; \; \; \; {\color[rgb]{0.064095,0.501831,0.501977}h} + \frac{1}{2}{\color[rgb]{0.059472,0.501943,0.998465}v}^2 = const \]
then
\[ \large \require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \]
and
\[ \large \require{color}\frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0. \]
Using the ideal gas law,
\[ \large {\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}, \]
the equation for enthalpy is
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}, \]
where
\[ \large \require{color}{\color[rgb]{0.986252,0.007236,0.027423}R}={\color[rgb]{0.501958,0.501942,0.014744}c_p}-c_v, \]
and
\[ \large \require{color}\gamma = \frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{c_v}. \]
From this, one obtains
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}= \frac{\gamma -1}{\gamma} {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.064095,0.501831,0.501977}h}. \]
Additionally, since we know that
\[ \large \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}} \]
and
\[ \large \require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a} \]
we can solve for
\[ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2, \]
yielding
\[ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.041893,0.355669,0.727621}M}^2{\color[rgb]{0.989013,0.435749,0.811750}a}^2={\color[rgb]{0.814433,0.253157,0.091125}\rho} \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 = {\color[rgb]{0.315209,0.728565,0.037706}p}\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2 \]
To expand
\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx} \]
we expand the derivative of
\[ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 \]
into terms of \({\color[rgb]{0.059472,0.501943,0.998465}v}\) and \({\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\) and we use the expression for \(p\) to expand the first term:
\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{\gamma -1}{\gamma}\left({\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.064095,0.501831,0.501977}h}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}\right) + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \] \[ \large = \frac{\gamma -1}{\gamma}\left({\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.064095,0.501831,0.501977}h}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}\right) + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \]
Using that \(\require{color}\frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \longrightarrow \frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} = -{\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}\) we can further expand the equation to get:
\[ \large \require{color} \frac{\gamma -1}{\gamma}\left(-{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + {\color[rgb]{0.064095,0.501831,0.501977}h}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}\right) + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \]
Then we can use that \(\require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}= \frac{\gamma -1}{\gamma} {\color[rgb]{0.064095,0.501831,0.501977}h}\) to further expand to
\[ \large \require{color}-\frac{\gamma -1}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \],
simplifying this we get:
\[ \large \require{color}\frac{1-\gamma}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + \frac{\gamma}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{1}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} \]
Then using the fact that \(\require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0\) we can get that
\[ \large \require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}={\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \] \[ \large \Rightarrow \frac{1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}=-\frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}, \]
subsituting this into the equation we get:
\[ \large \require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}={\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \longrightarrow \frac{1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}=-\frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \]
Lastly by using \(\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 = \frac{1}{\gamma}{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\) we can get that
\[ \large \require{color} \frac{1}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.059472,0.501943,0.998465}v}} \].
Substituting this in we get:
\[ \large \require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} - \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}({\color[rgb]{0.041893,0.355669,0.727621}M}^2 - 1)}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \]
Putting this back into the original force equation we have that:
\[ \large \require{color}-\frac{4{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w}{D} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}({\color[rgb]{0.041893,0.355669,0.727621}M}^2 - 1)}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \]
From dimensional analysis we know that
\[ \large \require{color}c_f = \frac{{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w}{\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}, \]
which upon substituting yields
\[ \large \require{color}\frac{4c_f}{D}dx = -\frac{2{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}({\color[rgb]{0.041893,0.355669,0.727621}M}^2-1)\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}}. \]
To solve the above, we need to integrate along \(x\) from \(0\) to the final length \(L\) and from the initial to final velocities:
\[ \large \require{color}\int_0^L \frac{4c_f}{D}dx = \int_{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}^{{\color[rgb]{0.059472,0.501943,0.998465}v}_2} -\frac{2{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}({\color[rgb]{0.041893,0.355669,0.727621}M}^2-1)\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} \]
By using that
\[ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 = \frac{}{\gamma}{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \]
we arrive at
\[ \large \require{color}-\frac{2{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2} = -\frac{2}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2}. \]
Substituting this in we get:
\[ \large \require{color}\int_0^L \frac{4c_f}{D}dx = \int_{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}^{{\color[rgb]{0.059472,0.501943,0.998465}v}_2} -\frac{2}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2}({\color[rgb]{0.041893,0.355669,0.727621}M}^2-1)\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} \]
We wish to integrate with one variable, in this case \({\color[rgb]{0.041893,0.355669,0.727621}M}^2\) as we want the relationship between Mach and coefficient of friction, and \({\color[rgb]{0.041893,0.355669,0.727621}M}^2\) is easier to integrate over. To do this we start with the equation for stagnation enthalpy
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.064095,0.501831,0.501977}h}+\frac{1}{2}{\color[rgb]{0.059472,0.501943,0.998465}v}^2 \]
and use
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T} \],
yielding
\[ \large \require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a} \],
and
\[ \large \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}} \]
to get:
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T} + \frac{1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} \]
Then using that \(\require{color}{\color[rgb]{0.986252,0.007236,0.027423}R}={\color[rgb]{0.501958,0.501942,0.014744}c_p}-c_v\) we can simplify to get:
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T} + {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2) \]
Next we solve for temperature from
\[ \large \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}} \]
\[ \large \Rightarrow {\color[rgb]{0.121820,0.954406,0.966585}T}=\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^2}{{\color[rgb]{0.986252,0.007236,0.027423}R}\gamma}=\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}} \],
plugging this in we have:
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}}(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2). \]
We need to integrate in order to get a relationship between \(d{\color[rgb]{0.059472,0.501943,0.998465}v}\) and \(d{\color[rgb]{0.041893,0.355669,0.727621}M}^2\); to make this easier we take the logarithmic of both sides and use logarithmic rules in order to separate the right hand side:
\[ \large \require{color}\textrm{ln}({\color[rgb]{0.064095,0.501831,0.501977}h}_0) = \textrm{ln}\left({\color[rgb]{0.501958,0.501942,0.014744}c_p}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}}\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)\right) \]
\[ \large \Rightarrow \textrm{ln} \left(\frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{{\color[rgb]{0.986252,0.007236,0.027423}R}\gamma}\right) + \textrm{ln}\left({\color[rgb]{0.059472,0.501943,0.998465}v}^2\right) - \textrm{ln}\left({\color[rgb]{0.041893,0.355669,0.727621}M}^2\right) + \textrm{ln}\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right) \]
Then we can take the derivative of both sides:
\[ \large \require{color}0 = 0 + \frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}(2{\color[rgb]{0.059472,0.501943,0.998465}v}d{\color[rgb]{0.059472,0.501943,0.998465}v}) - \frac{1}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2}d{\color[rgb]{0.041893,0.355669,0.727621}M}^2 + \frac{1}{1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2} \frac{\gamma-1}{2}d{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \]
\[ \large \Rightarrow \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} - \frac{d{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)} \]
Solving for \(\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\) we get that:
\[ \large \require{color}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} = \frac{d{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{2{\color[rgb]{0.041893,0.355669,0.727621}M}^2\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)} \]
Then by substituting this into the integral we have:
\[ \large \int_0^L \frac{4c_f}{D}dx = -\int_{{\color[rgb]{0.041893,0.355669,0.727621}M}_1}^{{\color[rgb]{0.041893,0.355669,0.727621}M}_2} \frac{1}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2} (1-{\color[rgb]{0.041893,0.355669,0.727621}M}^2) \frac{d{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\left(1+\frac{(\gamma-1)}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)} \]
Integrating with \({\color[rgb]{0.041893,0.355669,0.727621}M}_2 =1\) we get that:
\[ \large \require{color}\frac{4c_{f}}{D}L = \frac{1 - {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2} + \frac{\gamma +1}{2\gamma} \ \textrm{ln}\left(\frac{\frac{\gamma +1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}{1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}\right) \]