Shock wave solvers
If you are looking for a rapid way to compute all the ratios across the shock wave, you have come to the right place.
Two solvers are implemented: 1. shockwave_solver: use this to solve for normal or oblique shock waves. 2. conical_shockwave_solver: use this to solve for conical shock wave.
In the following, a few examples are presented to show how to use them and to understand some of the limitations.
[1]:
from pygasflow import (
shockwave_solver,
conical_shockwave_solver
)
shockwave_solver
Let’s look at the documentation for this function to see how it works:
[2]:
help(shockwave_solver)
Help on function shockwave_solver in module pygasflow.solvers.shockwave:
shockwave_solver(p1_name, p1_value, p2_name='beta', p2_value=90, gamma=1.4, flag='weak')
Try to compute all the ratios, angles and mach numbers across the shock wave.
Remember: a normal shock wave has a wave angle beta=90 deg.
Parameters
----------
p1_name : string
Name of the first parameter given in input. Can be either one of:
* ``'pressure'``: Pressure Ratio P2/P1
* ``'temperature'``: Temperature Ratio T2/T1
* ``'density'``: Density Ratio rho2/rho1
* ``'total_pressure'``: Total Pressure Ratio P02/P01
* ``'m1'``: Mach upstream of the shock wave
* ``'mn1'``: Normal Mach upstream of the shock wave
* ``'mn2'``: Normal Mach downstream of the shock wave
* ``'beta'``: The shock wave angle [in degrees]. It can only be used
if ``p2_name='theta'``.
* ``'theta'``: The deflection angle [in degrees]. It can only be
used if ``p2_name='beta'``.
If the parameter is a ratio, it is in the form downstream/upstream:
p1_value : float
Actual value of the parameter.
p2_name : string, optional
Name of the second parameter. It could either be:
* ``'beta'``: Shock wave angle.
* ``'theta'``: Flow deflection angle.
* ``'mn1'``: Input Normal Mach number.
Default to ``'beta'``.
p2_value : float, optional
Value of the angle in degrees.
Default to 90 degrees (normal shock wave).
gamma : float, optional
Specific heats ratio. Default to 1.4. Must be > 1.
flag : string, optional
Chose what solution to compute if the angle 'theta' is provided.
Can be either ``'weak'`` or ``'strong'``. Default to ``'weak'``.
Returns
-------
M1 : float
Mach number upstream of the shock wave.
Mn1 : float
Normal Mach number upstream of the shock wave.
M2 : float
Mach number downstream of the shock wave.
Mn2 : float
Normal Mach number downstream of the shock wave.
beta : float
Shock wave angle in degrees.
theta : float
Flow deflection angle in degrees.
pr : float
Pressure ratio across the shock wave.
dr : float
Density ratio across the shock wave.
tr : float
Temperature ratio across the shock wave.
tpr : float
Total Pressure ratio across the shock wave.
As we can see, it is really quite simple to use it. Just insert the name of the known parameter, its value, the name of the angle and its value, the specific heats ratio and you are good to go.
Let’s first create a simple function to print the results:
[3]:
def print_oblique_shockwave(M1, MN1, M2, MN2, beta, theta, pr, dr, tr, tpr):
print("M1 \t\t {}".format(M1))
print("Mn1 \t\t {}".format(MN1))
print("M2 \t\t {}".format(M2))
print("Mn2 \t\t {}".format(MN2))
print("beta \t\t {}".format(beta))
print("theta \t\t {}".format(theta))
print("p2/p1 \t\t {}".format(pr))
print("rho2/rho1 \t {}".format(dr))
print("t2/t1 \t\t {}".format(tr))
print("p02/p01 \t {}".format(tpr))
print()
Normal Shock Wave
By default, this function can be used to solve normal shock wave on air (specific heats ratio \(\gamma=1.4\) and shock wave angle \(\beta=90\)):
[4]:
result = shockwave_solver("m1", 4)
print_oblique_shockwave(*result)
M1 4.0
Mn1 4.0
M2 0.43495883620084
Mn2 0.43495883620084
beta 90.0
theta 1.2529838033097996e-14
p2/p1 18.5
rho2/rho1 4.571428571428572
t2/t1 4.046874999999999
p02/p01 0.13875622089168838
Because this function is designed to solve also oblique shock waves, we have some parameters that we can discard: 1. Mn1, the normal Mach number upstream of the shock wave, which is obviously equal to M1. 2. Mn2, the normal Mach number downstream of the shock wave, which is obviously equal to M2. 3. beta, the shock wave angle, obviously \(90°\) since we are solving a normal shock wave. 4. theta, the deflection angle. We notice a very small value. Indeed, in this case it should be 0, but numerical errors are preventing this result. Anyway, since we are solving a normal shock wave, we do not need to worry about this angle.
Let’s look at the same example, but inserting different parameters. We expect the very same results, and indeed we get them.
[5]:
r1 = shockwave_solver("pressure", 18.5)
r2 = shockwave_solver("density", 4.571428571428572)
r3 = shockwave_solver("temperature", 4.046874999999999)
r4 = shockwave_solver("mn2", 0.43495883620084)
print_oblique_shockwave(*r1)
print_oblique_shockwave(*r2)
print_oblique_shockwave(*r3)
print_oblique_shockwave(*r4)
M1 4.0
Mn1 4.0
M2 0.43495883620084
Mn2 0.43495883620084
beta 90
theta 1.2529838033097996e-14
p2/p1 18.5
rho2/rho1 4.571428571428572
t2/t1 4.046874999999999
p02/p01 0.13875622089168838
M1 4.0
Mn1 4.0
M2 0.43495883620084
Mn2 0.43495883620084
beta 90
theta 1.2529838033097996e-14
p2/p1 18.5
rho2/rho1 4.571428571428572
t2/t1 4.046874999999999
p02/p01 0.13875622089168838
M1 4.0
Mn1 4.0
M2 0.43495883620084
Mn2 0.43495883620084
beta 90
theta 1.2529838033097996e-14
p2/p1 18.5
rho2/rho1 4.571428571428572
t2/t1 4.046874999999999
p02/p01 0.13875622089168838
M1 4.0
Mn1 4.0
M2 0.43495883620084
Mn2 0.43495883620084
beta 90
theta 1.2529838033097996e-14
p2/p1 18.5
rho2/rho1 4.571428571428572
t2/t1 4.046874999999999
p02/p01 0.13875622089168838
What if I want to solve a normal shock wave on a gas different than air? Obviously I need to change the specific heats ratio:
[6]:
result = shockwave_solver("m1", 4, gamma=1.2)
print_oblique_shockwave(*result)
M1 4.0
Mn1 4.0
M2 0.3689521031926255
Mn2 0.3689521031926255
beta 90.0
theta 2.024050759192753e-14
p2/p1 17.36363636363636
rho2/rho1 6.769230769230771
t2/t1 2.5650826446280983
p02/p01 0.06095825251291948
Oblique Shock Wave
Let’s now consider an oblique shock wave. Let’s start with air, upstream Mach number \(M_{1} = 4\) and a deflection angle \(\theta = 20°\).
[7]:
result = shockwave_solver("m1", 4, 'theta', 20)
print_oblique_shockwave(*result)
M1 4.0
Mn1 2.147072259523833
M2 2.5686168890322807
Mn2 0.5543701693902555
beta 32.46389685027039
theta 20.0
p2/p1 5.211572502219574
rho2/rho1 2.8782256018884964
t2/t1 1.8106893701453053
p02/p01 0.6524015014542756
As we have seen from the documentation, if the deflection angle \(\theta\) is provided, by default the function will solve for the *weak shock solution*. If we want to solve for the *strong shock solution*, we must change the flag argument. But in doing so, we also need to specify the specific heats ratio.
[8]:
result = shockwave_solver("m1", 4, 'theta', 20, gamma=1.4, flag='strong')
print_oblique_shockwave(*result)
M1 4.0
Mn1 3.977010652739413
M2 0.48523272499938425
Mn2 0.43558154682352673
beta 83.85418932598304
theta 20.0
p2/p1 18.286049354003236
rho2/rho1 4.5588434129476605
t2/t1 4.011115912002739
p02/p01 0.14147889026890678
Note that the solver is able to detect detachment. Let’s say we provide a flow deflection angle greater that the maximum allowed for the provided upstream Mach number (see Oblique Shock Diagram):
[9]:
result = shockwave_solver("m1", 2, 'theta', 30)
print_oblique_shockwave(*result)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [9], in <module>
----> 1 result = shockwave_solver("m1", 2, 'theta', 30)
2 print_oblique_shockwave(*result)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:206, in as_array.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
204 if i < len(args):
205 args[i] = convert_to_ndarray(args[i])
--> 206 return original_function(*args, **kwargs)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:121, in check_shockwave.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
118 if "flag" in all_param.keys() and len(args) == len(all_param.keys()):
119 args[-1] = _check_flag_shockwave(all_param["flag"])
--> 121 res = original_function(*args, **kwargs)
122 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/solvers/shockwave.py:163, in shockwave_solver(p1_name, p1_value, p2_name, p2_value, gamma, flag)
160 M2 = MN2 / np.sin(np.deg2rad(beta - theta))
161 elif M1 is not None:
162 # at this point, either beta or theta is set, not both!
--> 163 MN1 = normal_mach_upstream.__no_check(M1, beta, theta, gamma, flag)
164 # compute the different ratios
165 pr, dr, tr, tpr, MN2 = get_ratios_from_normal_mach_upstream.__no_check(MN1, gamma)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:125, in check_shockwave.<locals>.decorator.<locals>.no_check_function(*args, **kwargs)
124 def no_check_function(*args, **kwargs):
--> 125 res = original_function(*args, **kwargs)
126 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/shockwave.py:497, in normal_mach_upstream(M1, beta, theta, gamma, flag)
495 theta_max = max_theta_from_mach(M1, gamma)
496 if np.any(theta > theta_max):
--> 497 raise ValueError("Detachment detected: can't solve the flow when theta > theta_max.\n" +
498 "M1 = {}\n".format(M1) +
499 "theta_max(M1) = {}\n".format(theta_max) +
500 "theta = {}\n".format(theta))
501 beta = beta_from_mach_theta(M1, theta)
502 MN1 = dict()
ValueError: Detachment detected: can't solve the flow when theta > theta_max.
M1 = [2.]
theta_max(M1) = 22.97353176093536
theta = [30.]
[10]:
print_oblique_shockwave(*shockwave_solver("m1", 1.8400645259024768, "theta", 20))
M1 1.8400645259024768
Mn1 1.664726146515799
M2 0.9225790066302748
Mn2 0.64990193322869
beta 64.78435329607117
theta 20.0
p2/p1 3.066532000042231
rho2/rho1 2.1396485447978377
t2/t1 1.433194254027345
p02/p01 0.8701212937216544
Again, we can insert different parameter to solve for the shock. For example, let’s say I know the pressure ratio across the shock, the deflection angle, and I want to solve for the *weak shock solution*. I can do:
[11]:
result = shockwave_solver("pressure", 5.211572502219574, 'theta', 20)
print_oblique_shockwave(*result)
M1 nan
Mn1 2.147072259523833
M2 nan
Mn2 0.5543701693902555
beta nan
theta 20.0
p2/p1 5.211572502219574
rho2/rho1 2.8782256018884964
t2/t1 1.8106893701453053
p02/p01 0.6524015014542756
/home/davide/Documents/Development/pygasflow/pygasflow/solvers/shockwave.py:194: UserWarning: Undetermined case. Setting M1 = beta = M2 = NaN
warnings.warn("Undetermined case. Setting M1 = beta = M2 = NaN")
Comparing with the previous results, we see that the ratios are correct, the normal Mach numbers are correct, but the shock wave angle \(\beta\) and thus the Mach numbers upstream and downstream of the shock are not determined. This is to be expected: given a ratio, we can compute only the normal Mach number across the shock wave corresponsing to the weak solution, but we are unable to uniquely determine the shock wave angle \(\beta\) and the Mach number.
Therefore, if we would like to compute the *strong shock solution* giving a ratio and the deflection angle :math:`theta`, we will get the wrong results:
[12]:
result = shockwave_solver("pressure", 5.211572502219574, 'theta', 20, gamma=1.4, flag='strong')
print_oblique_shockwave(*result)
M1 nan
Mn1 2.147072259523833
M2 nan
Mn2 0.5543701693902555
beta nan
theta 20.0
p2/p1 5.211572502219574
rho2/rho1 2.8782256018884964
t2/t1 1.8106893701453053
p02/p01 0.6524015014542756
Even by specify the flag='strong', at the moment of writing this tutorial, it is not possible to compute the Normal Mach number corresponding to strong solution of the input ration, and it is not possible to uniquely determine the shock wave angle \(\beta\) and thus computing the correct results.
We can solve a shock wave by knowing both the shock wave angle \(\beta\) and the deflection angle \(\theta\). Following the previus example, by chosing the shock wave angle corresponding to the weak solution we get:
[14]:
r1 = shockwave_solver("beta", 32.46389685027039, 'theta', 20)
r2 = shockwave_solver('theta', 20, "beta", 32.46389685027039)
print_oblique_shockwave(*r1)
print()
print_oblique_shockwave(*r2)
M1 4.000000000000003
Mn1 2.1470722595238345
M2 2.5686168890322807
Mn2 0.5543701693902553
beta 32.46389685027039
theta 20.000000000000004
p2/p1 5.211572502219582
rho2/rho1 2.8782256018884977
t2/t1 1.8106893701453068
p02/p01 0.6524015014542744
M1 4.000000000000003
Mn1 2.1470722595238345
M2 2.5686168890322807
Mn2 0.5543701693902553
beta 32.46389685027039
theta 20.000000000000004
p2/p1 5.211572502219582
rho2/rho1 2.8782256018884977
t2/t1 1.8106893701453068
p02/p01 0.6524015014542744
Note that we do not need to specify the flag arguments, since for each \((\beta, \theta)\) corresponds only one Mach number.
By chosing the shock wave angle corresponding to the strong solution we get:
[15]:
r1 = shockwave_solver("beta", 83.85418932598304, 'theta', 20)
r2 = shockwave_solver('theta', 20, "beta", 83.85418932598304)
print_oblique_shockwave(*r1)
print_oblique_shockwave(*r2)
M1 3.9999999999999862
Mn1 3.9770106527393994
M2 0.48523272499938463
Mn2 0.4355815468235271
beta 83.85418932598304
theta 19.99999999999999
p2/p1 18.286049354003108
rho2/rho1 4.558843412947653
t2/t1 4.011115912002717
p02/p01 0.14147889026890847
M1 3.9999999999999862
Mn1 3.9770106527393994
M2 0.48523272499938463
Mn2 0.4355815468235271
beta 83.85418932598304
theta 19.99999999999999
p2/p1 18.286049354003108
rho2/rho1 4.558843412947653
t2/t1 4.011115912002717
p02/p01 0.14147889026890847
Finally, the solver should be able to determine if a solution doesn’t exist for the current choice of parameters. For example (see Oblique Shock Diagram):
[16]:
print_oblique_shockwave(*shockwave_solver("beta", 88, "theta", 20))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [16], in <module>
----> 1 print_oblique_shockwave(*shockwave_solver("beta", 88, "theta", 20))
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:206, in as_array.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
204 if i < len(args):
205 args[i] = convert_to_ndarray(args[i])
--> 206 return original_function(*args, **kwargs)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:121, in check_shockwave.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
118 if "flag" in all_param.keys() and len(args) == len(all_param.keys()):
119 args[-1] = _check_flag_shockwave(all_param["flag"])
--> 121 res = original_function(*args, **kwargs)
122 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/solvers/shockwave.py:150, in shockwave_solver(p1_name, p1_value, p2_name, p2_value, gamma, flag)
148 if not isinstance(theta, np.ndarray):
149 theta = theta * np.ones_like(beta)
--> 150 M1 = mach_from_theta_beta.__no_check(theta, beta)
151 else: # 'm2'
152 # TODO:
153 # Is it even possible to solve it knowing only M2, beta or M2, theta?????
154 raise NotImplementedError("Solving a shock wave with a given M2 is not yet implemented.")
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:125, in check_shockwave.<locals>.decorator.<locals>.no_check_function(*args, **kwargs)
124 def no_check_function(*args, **kwargs):
--> 125 res = original_function(*args, **kwargs)
126 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/shockwave.py:800, in mach_from_theta_beta(theta, beta, gamma)
798 theta_max[i] = theta_from_mach_beta(100000, b, gamma)
799 if np.any(theta > theta_max):
--> 800 raise ValueError("There is no solution for the current choice of" +
801 " parameters. Please check the Oblique Shock diagram with the following"
802 " parameters:\n" +
803 "beta = {}\n".format(beta) +
804 "theta = {}\n".format(theta))
805 # case beta == 90 and theta == 0, which leaves M to be indeterminate, NaN
806 idx0 = np.bitwise_and(beta == 90, theta == 0)
ValueError: There is no solution for the current choice of parameters. Please check the Oblique Shock diagram with the following parameters:
beta = [88.]
theta = [20.]
conical_shockwave_solver
Let’s look at the documentation for this function to see how it works:
[17]:
help(conical_shockwave_solver)
Help on function conical_shockwave_solver in module pygasflow.solvers.shockwave:
conical_shockwave_solver(M1, param_name, param_value, gamma=1.4, flag='weak')
Try to compute all the ratios, angles and mach numbers across the conical shock wave.
Parameters
----------
M1 : float
Upstream Mach number. Must be M1 > 1.
param_name : string
Name of the parameter given in input. Can be either one of:
* ``'mc'``: Mach number at the cone's surface.
* ``'theta_c'``: Half cone angle.
* ``'beta'``: shock wave angle.
param_value : float
Actual value of the parameter. Requirements:
* Mc >= 0
* 0 < beta <= 90
* 0 < theta_c < 90
gamma : float, optional
Specific heats ratio. Default to 1.4. Must be > 1.
flag : string, optional
Can be either ``'weak'`` or ``'strong'``. Default to ``'weak'``
(in conical shockwaves, the strong solution is rarely encountered).
Returns
-------
M : float
Upstream Mach number.
Mc : float
Mach number at the surface of the cone.
theta_c : float
Half cone angle.
beta : float
Shock wave angle.
delta : float
Flow deflection angle.
pr : float
Pressure ratio across the shock wave.
dr : float
Density ratio across the shock wave.
tr : float
Temperature ratio across the shock wave.
tpr : float
Total Pressure ratio across the shock wave.
pc_p1 : float
Pressure ratio between the cone's surface and the upstream condition.
rhoc_rho1 : float
Density ratio between the cone's surface and the upstream condition.
Tc_T1 : float
Temperature ratio between the cone's surface and the upstream condition.
[18]:
def print_conical_shockwave(M1, Mc, theta_c, beta, delta, pr, dr, tr, tpr, pc_p1, rhoc_rho1, Tc_T1):
print("M1 \t\t {}".format(M1))
print("Mc \t\t {}".format(Mc))
print("theta_c \t {}".format(theta_c))
print("beta \t\t {}".format(beta))
print("delta \t\t {}".format(delta))
print("p2/p1 \t\t {}".format(pr))
print("rho2/rho1 \t {}".format(dr))
print("t2/t1 \t\t {}".format(tr))
print("p02/p01 \t {}".format(tpr))
print("pc/p1 \t\t {}".format(pc_p1))
print("rho_c/rho1 \t {}".format(rhoc_rho1))
print("tc/t1 \t\t {}".format(Tc_T1))
[19]:
result = conical_shockwave_solver(4, 'beta', 20)
print_conical_shockwave(*result)
M1 4.0
Mc 3.3494441585129917
theta_c 12.923342940926778
beta 20
delta 7.444220294251439
p2/p1 2.0169185308895377
rho2/rho1 1.6342327959123613
t2/t1 1.2341684342245318
p02/p01 0.965779172840477
pc/p1 2.385514639007803
rho_c/rho1 1.8423870467954182
tc/t1 1.2947955985454198
[20]:
result = conical_shockwave_solver(3, 'theta_c', 26.377485151514648)
print_conical_shockwave(*result)
M1 3.0
Mc 2.000000000000039
theta_c 26.377485151514648
beta 35.90847523748602
delta 18.392639366690542
p2/p1 3.44505054944179
rho2/rho1 2.294355459847854
t2/t1 1.5015330491423686
p02/p01 0.8304690737856829
pc/p1 3.898726823356718
rho_c/rho1 2.506324386443648
tc/t1 1.5555555555555283
[21]:
result = conical_shockwave_solver(3, 'mc', 2)
print_conical_shockwave(*result)
M1 3.0
Mc 2.0
theta_c 26.377485151514648
beta 35.90847523748602
delta 18.392639366690542
p2/p1 3.44505054944179
rho2/rho1 2.294355459847854
t2/t1 1.5015330491423686
p02/p01 0.8304690737856829
pc/p1 3.8987268233569545
rho_c/rho1 2.5063243864437568
tc/t1 1.5555555555555554
[22]:
result = conical_shockwave_solver(3, 'beta', 35.90847523748602)
print_conical_shockwave(*result)
M1 3.0
Mc 2.000000000000039
theta_c 26.377485151514648
beta 35.90847523748602
delta 18.392639366690542
p2/p1 3.44505054944179
rho2/rho1 2.294355459847854
t2/t1 1.5015330491423686
p02/p01 0.8304690737856829
pc/p1 3.898726823356718
rho_c/rho1 2.506324386443648
tc/t1 1.5555555555555283
The solver should be able to detect detachment:
[23]:
result = conical_shockwave_solver(2, 'theta_c', 45)
print_conical_shockwave(*result)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [23], in <module>
----> 1 result = conical_shockwave_solver(2, 'theta_c', 45)
2 print_conical_shockwave(*result)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:206, in as_array.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
204 if i < len(args):
205 args[i] = convert_to_ndarray(args[i])
--> 206 return original_function(*args, **kwargs)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:121, in check_shockwave.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
118 if "flag" in all_param.keys() and len(args) == len(all_param.keys()):
119 args[-1] = _check_flag_shockwave(all_param["flag"])
--> 121 res = original_function(*args, **kwargs)
122 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/solvers/shockwave.py:285, in conical_shockwave_solver(M1, param_name, param_value, gamma, flag)
283 Mc, theta_c = mach_cone_angle_from_shock_angle(M1, beta, gamma)
284 elif theta_c:
--> 285 Mc, _, beta = shock_angle_from_mach_cone_angle(M1, theta_c, gamma, flag)
287 # compute the ratios across the shockwave
288 MN1 = M1 * np.sin(np.deg2rad(beta))
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:206, in as_array.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
204 if i < len(args):
205 args[i] = convert_to_ndarray(args[i])
--> 206 return original_function(*args, **kwargs)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:121, in check_shockwave.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
118 if "flag" in all_param.keys() and len(args) == len(all_param.keys()):
119 args[-1] = _check_flag_shockwave(all_param["flag"])
--> 121 res = original_function(*args, **kwargs)
122 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/shockwave.py:1161, in shock_angle_from_mach_cone_angle(M1, theta_c, gamma, flag)
1159 Mc = np.zeros_like(M1)
1160 for i, m in enumerate(M1):
-> 1161 Mc[i], theta_c_comp[i], beta[i] = function(m)
1162 return ret_correct_vals(Mc), ret_correct_vals(theta_c_comp), ret_correct_vals(beta)
1163 return function(M1)
File ~/Documents/Development/pygasflow/pygasflow/shockwave.py:1135, in shock_angle_from_mach_cone_angle.<locals>.function(M)
1133 Mc, tcmax, bmax = max_theta_c_from_mach(M)
1134 if theta_c > tcmax:
-> 1135 raise ValueError("Detachment detected: can't solve the flow when theta_c > theta_c_max.\n" +
1136 "M1 = {}\n".format(M1) +
1137 "theta_c_max(M1) = {}\n".format(tcmax) +
1138 "theta_c = {}\n".format(theta_c))
1140 # need to add a small offset otherwise mach_cone_angle_from_shock_angle
1141 # could crash
1142 bmin = np.rad2deg(np.arcsin(1 / M)) + 1e-08
ValueError: Detachment detected: can't solve the flow when theta_c > theta_c_max.
M1 = [2.]
theta_c_max(M1) = 40.688476890932144
theta_c = 45
[ ]: