```python
import numpy as np
from scipy import optimize
```



Starting with an initial flow over a horizontal surface where $M_1 = 2.2$, an oblique shock forms at an angle of $\theta = 35^{\circ}$, which deflects the flow by the angle $\delta$. The flow now has Mach number $M_2$.

To satisfy the boundary condition of the surface, the oblique shock reflects at an angle $\beta$ from the surface, which deflects the flow back to the horizontal direction with a flow at a Mach number of $M_3$.

What is the angle of reflection ($\beta$)? How do the strengths of each shock compare?

First, let's determine the deflection angle of the first oblique shock using
\begin{equation}
\tan \delta = 2 (\cot \theta) \left[ \frac{M_1^2 \sin^2 \theta - 1}{M_1^2 (\gamma + \cos 2\theta) + 2}\right]
\end{equation}


```python
# known values
gamma = 1.4
M1 = 2.2
theta = np.radians(35.0) # convert from degrees to radians

delta = np.arctan(
    2 * (1 / np.tan(theta)) * (M1**2 * np.sin(theta)**2 - 1) / (M1**2 * (gamma + np.cos(2*theta)) + 2)
)
print(f'Deflection angle: {np.degrees(delta):.3}')
```

    Deflection angle: 9.21


Now we can determine the conditions after the oblique shock, using
\begin{equation}
M_{1n} = M_1 \sin \theta
\end{equation}
and
\begin{equation}
M_{2n}^2 = \frac{M_{1n}^2 + \frac{2}{\delta-1}}{\frac{2\gamma}{\gamma-1} M_{1n}^2 - 1}
\end{equation}


```python
M1n = M1 * np.sin(theta)
M2n = np.sqrt((M1n**2 + (2/(gamma - 1))) / (M1n**2 * (2*gamma)/(gamma - 1) - 1))

M2 = M2n / np.sin(theta - delta)
print(f'M_2 = {M2:.3}')
```

    M_2 = 1.85


Now, the reflected oblique shock must turn the flow back to the horizontal direction to satisfy flow tangency with the wall, so $\delta_2 = \delta$. We can thus determine the angle of the second shock, $\theta_2$:


```python
def oblique(theta, M1, delta, gamma):
    return (
        np.tan(delta) - 
        2 * (1 / np.tan(theta)) * (M1**2 * np.sin(theta)**2 - 1) / (M1**2 * (gamma + np.cos(2*theta)) + 2)
    )

root = optimize.root_scalar(oblique, # function we are solving
                            args=(M2, delta, gamma),
                            # give range for weak shocks
                            bracket=[np.radians(0.0001), np.radians(45.0)]
                           )
theta_2 = root.root
print(f'theta_2 = {np.degrees(theta_2):.3} degrees')
```

    theta_2 = 41.7 degrees


 **But** this shock angle is defined with respect to the flow direction in the middle region, so $\beta = \theta_2 - \delta$:


```python
beta = theta_2 - delta
print(f'beta = {np.degrees(beta):.3} degrees')
```

    beta = 32.5 degrees


Thus, the angle of reflection $\beta$ is **smaller** than the angle of incidence $\theta = 35^{\circ}$. Regarding the shock strength, we can compare the normal Mach numbers of the flow prior to each shock:


```python
print(f'M1n = {M1n:.4}')
M3n = M2 * np.sin(theta_2)
print(f'M2n = {M3n:.4}')
```

    M1n = 1.262
    M2n = 1.233


We can see that the second shock occurs at a smaller Mach number, so it is **weaker**.

f43ce2b9dded71a17c27009fd4a2c533689b90f3
ipynb
Jupyter Notebook
oblique_shock.ipynb
kyleniemeyer/gasdynamics
50dce7a030daa2757aa55cf6c9a5ae66d079ab72
MIT
2019-10-17T18:21:23.000Z
2021-08-17T19:30:07.000Z
oblique_shock.ipynb
kyleniemeyer/gasdynamics
50dce7a030daa2757aa55cf6c9a5ae66d079ab72
MIT
oblique_shock.ipynb
kyleniemeyer/gasdynamics
50dce7a030daa2757aa55cf6c9a5ae66d079ab72
MIT
Qwen/Qwen-72B
1. YES
2. YES