Fanno Solver
Fanno ratios plot
Before explaining the solver, just a quick visual sanity check with the ratios plot.
[1]:
import numpy as np
import matplotlib.pyplot as plt
from pygasflow import fanno_solver
gamma = 1.4
M = np.linspace(1e-05, 5, 100)
results = fanno_solver("m", M, gamma)
plt.figure()
plt.plot(M, results[1], label=r"$p / p^{*}$")
plt.plot(M, results[2], label=r"$\rho / \rho^{*}$")
plt.plot(M, results[3], label=r"$T / T^{*}$")
plt.plot(M, results[4], label=r"$P0 / P_{0}^{*}$")
plt.plot(M, results[5], label=r"$U / U^{*}$")
plt.plot(M, results[6], label=r"$4f L^{*} / D$")
plt.plot(M, results[7], label=r"$(s^{*}-s)/R$")
plt.xlim(0, max(M))
plt.ylim(0, 3)
plt.xlabel("M")
plt.ylabel("ratios")
plt.grid(which='major', linestyle='-', alpha=0.7)
plt.grid(which='minor', linestyle=':', alpha=0.5)
plt.legend(loc='upper right')
plt.suptitle("Fanno Flow")
plt.title("$\gamma$ = " + str(gamma))
plt.show()
Examples
This solver allows to quickly compute the Mach number, all the ratios and parameters (entropy paramater, critical friction) by specifying what parameter we do know.
Let’s read the solver’s documentation:
[2]:
help(fanno_solver)
Help on function fanno_solver in module pygasflow.solvers.fanno:
fanno_solver(param_name, param_value, gamma=1.4, to_dict=False)
Compute all Fanno ratios and Mach number given an input parameter.
Parameters
----------
param_name : string
Name of the parameter given in input. Can be either one of:
* ``'m'``: Mach number
* ``'pressure'``: Critical Pressure Ratio P/P*
* ``'density'``: Critical Density Ratio rho/rho*
* ``'temperature'``: Critical Temperature Ratio T/T*
* ``'total_pressure_sub'``: Critical Total Pressure Ratio P0/P0*
for subsonic case.
* ``'total_pressure_super'``: Critical Total Pressure Ratio P0/P0*
for supersonic case.
* ``'velocity'``: Critical Velocity Ratio U/U*.
* ``'friction_sub'``: Critical Friction parameter 4fL*/D for
subsonic case.
* ``'friction_super'``: Critical Friction parameter 4fL*/D for
supersonic case.
* ``'entropy_sub'``: Entropy parameter (s*-s)/R for subsonic case.
* ``'entropy_super'``: Entropy parameter (s*-s)/R for supersonic case.
param_value : float/list/array_like
Actual value of the parameter. If float, list, tuple is given as
input, a conversion will be attempted.
gamma : float, optional
Specific heats ratio. Default to 1.4. Must be > 1.
to_dict : bool, optional
If False, the function returns a list of results. If True, it returns
a dictionary in which the keys are listed in the Returns section.
Default to False (return a list of results).
Returns
-------
m : array_like
Mach number
prs : array_like
Critical Pressure Ratio P/P*
drs : array_like
Critical Density Ratio rho/rho*
trs : array_like
Critical Temperature Ratio T/T*
tprs : array_like
Critical Total Pressure Ratio P0/P0*
urs : array_like
Critical Velocity Ratio U/U*
fps : array_like
Critical Friction Parameter 4fL*/D
eps : array_like
Critical Entropy Ratio (s*-s)/R
Examples
--------
Compute all ratios starting from a Mach number:
>>> from pygasflow import fanno_solver
>>> fanno_solver("m", 2)
[2.0, 0.408248290463863, 0.6123724356957945, 0.6666666666666667, 1.6875000000000002, 1.632993161855452, 0.3049965025814798, 0.523248143764548]
Compute the subsonic Mach number starting from the critical friction
parameter:
>>> results = fanno_solver("friction_sub", 0.3049965025814798)
>>> print(results[0])
0.6572579935727846
Compute the critical temperature ratio starting from multiple Mach numbers
for a gas having specific heat ratio gamma=1.2:
>>> results = fanno_solver("m", [0.5, 1.5], 1.2)
>>> print(results[3])
[1.07317073 0.89795918]
Compute the critical temperature ratio starting from multiple Mach numbers
for a gas having specific heat ratio gamma=1.2, returning a dictionary:
>>> results = fanno_solver("m", [0.5, 1.5], 1.2, to_dict=True)
>>> print(results["trs"])
[1.07317073 0.89795918]
This is just a pretty print function…
[3]:
def print_fanno(M, prs, drs, trs, tprs, urs, fps, eps):
print("M \t\t {}".format(M))
print("P/P* \t\t {}".format(prs))
print("rho/rho* \t {}".format(drs))
print("T/T* \t\t {}".format(trs))
print("P0/P0* \t\t {}".format(tprs))
print("U/U* \t\t {}".format(urs))
print("4fL*/D \t\t {}".format(fps))
print("(s*-s)/R \t {}".format(eps))
print()
[4]:
result = fanno_solver('m', 2)
print_fanno(*result)
M 2.0
P/P* 0.408248290463863
rho/rho* 0.6123724356957945
T/T* 0.6666666666666667
P0/P0* 1.6875000000000002
U/U* 1.632993161855452
4fL*/D 0.3049965025814798
(s*-s)/R 0.523248143764548
[5]:
r1 = fanno_solver('pressure', 0.408248290463863)
r2 = fanno_solver('density', 0.6123724356957945)
r3 = fanno_solver('temperature', 0.6666666666666667)
r4 = fanno_solver('total_pressure_super', 1.6875000000000002)
r5 = fanno_solver('velocity', 1.632993161855452)
r6 = fanno_solver('friction_super', 0.3049965025814798)
r7 = fanno_solver('entropy_super', 0.523248143764548)
print_fanno(*r1)
print_fanno(*r2)
print_fanno(*r3)
print_fanno(*r4)
print_fanno(*r5)
print_fanno(*r6)
print_fanno(*r7)
M 2.0000000000000004
P/P* 0.40824829046386285
rho/rho* 0.6123724356957945
T/T* 0.6666666666666665
P0/P0* 1.6875000000000007
U/U* 1.6329931618554523
4fL*/D 0.3049965025814797
(s*-s)/R 0.5232481437645482
M 2.0
P/P* 0.408248290463863
rho/rho* 0.6123724356957945
T/T* 0.6666666666666667
P0/P0* 1.6875000000000002
U/U* 1.632993161855452
4fL*/D 0.3049965025814798
(s*-s)/R 0.523248143764548
M 2.0000000000000004
P/P* 0.40824829046386285
rho/rho* 0.6123724356957945
T/T* 0.6666666666666665
P0/P0* 1.6875000000000007
U/U* 1.6329931618554523
4fL*/D 0.3049965025814797
(s*-s)/R 0.5232481437645482
M 1.9999999999986215
P/P* 0.40824829046426947
rho/rho* 0.6123724356960291
T/T* 0.6666666666670751
P0/P0* 1.6874999999980616
U/U* 1.6329931618548268
4fL*/D 0.30499650258106925
(s*-s)/R 0.5232481437633991
M 2.0
P/P* 0.408248290463863
rho/rho* 0.6123724356957945
T/T* 0.6666666666666667
P0/P0* 1.6875000000000002
U/U* 1.632993161855452
4fL*/D 0.3049965025814798
(s*-s)/R 0.523248143764548
M 1.9999999999986215
P/P* 0.40824829046426947
rho/rho* 0.6123724356960291
T/T* 0.6666666666670751
P0/P0* 1.6874999999980616
U/U* 1.6329931618548268
4fL*/D 0.30499650258106925
(s*-s)/R 0.5232481437633991
M 1.9999999999986215
P/P* 0.40824829046426947
rho/rho* 0.6123724356960291
T/T* 0.6666666666670751
P0/P0* 1.6874999999980616
U/U* 1.6329931618548268
4fL*/D 0.30499650258106925
(s*-s)/R 0.5232481437633991
Should you wish to solve for more than one Mach number of ratio at the same time, you can do:
[6]:
result = fanno_solver('m', [0.5, 1, 2, 4])
print_fanno(*result)
M [0.5 1. 2. 4. ]
P/P* [2.13808994 1. 0.40824829 0.13363062]
rho/rho* [1.87082869 1. 0.61237244 0.46770717]
T/T* [1.14285714 1. 0.66666667 0.28571429]
P0/P0* [ 1.33984375 1. 1.6875 10.71875 ]
U/U* [0.53452248 1. 1.63299316 2.13808994]
4fL*/D [1.06906031 0. 0.3049965 0.63306493]
(s*-s)/R [0.292553 0. 0.52324814 2.37199454]
We can also use this approach for other parameters. Like this:
[7]:
result = fanno_solver('friction_super', [1.06906031, 0, 0.3049965, 0.63306493])
print_fanno(*result)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [7], in <module>
----> 1 result = fanno_solver('friction_super', [1.06906031, 0, 0.3049965, 0.63306493])
2 print_fanno(*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:171, in check.<locals>.decorator.<locals>.wrapper_function(*args, **kwargs)
169 else:
170 kwargs["flag"] = flag
--> 171 res = original_function(*args, **kwargs)
172 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/solvers/fanno.py:115, in fanno_solver(param_name, param_value, gamma, to_dict)
113 M = fanno.m_from_critical_friction.__no_check(param_value, "sub", gamma)
114 elif param_name == 'friction_super':
--> 115 M = fanno.m_from_critical_friction.__no_check(param_value, "super", gamma)
116 elif param_name == 'entropy_sub':
117 M = fanno.m_from_critical_entropy.__no_check(param_value, "sub", gamma)
File ~/Documents/Development/pygasflow/pygasflow/utils/decorators.py:175, in check.<locals>.decorator.<locals>.no_check_function(*args, **kwargs)
174 def no_check_function(*args, **kwargs):
--> 175 res = original_function(*args, **kwargs)
176 return ret_correct_vals(res)
File ~/Documents/Development/pygasflow/pygasflow/fanno.py:329, in m_from_critical_friction(fp, flag, gamma)
327 upper_lim = ((gamma + 1) * np.log((gamma + 1) / (gamma - 1)) - 2) / (2 * gamma)
328 if np.any(fp < 0) or np.any(fp > upper_lim):
--> 329 raise ValueError('It must be 0 <= fp <= {}'.format(upper_lim))
331 # TODO: when solving the supersonic case, and ratio -> upper limit,
332 # I get: ValueError: f(a) and f(b) must have different signs
333 # need to be dealt with!
334
335 # func = lambda M, r: r - Critical_Friction_Parameter.__no_check(M, gamma)
336 func = lambda M, r: r - (((gamma + 1) / (2 * gamma)) * np.log(((gamma + 1) / 2) * M**2 / (1 + ((gamma - 1) / 2) * M**2)) + (1 / gamma) * (1 / M**2 - 1))
ValueError: It must be 0 <= fp <= 0.8215081164811902
Are you wondering what happened? The first two values inserted correspond to the subsonic solution, but we have asked the solver for the supersonic solution. At those values there is none. Therefore, attention must be used with all the parameteres that use different formulations for subsonic and supersonic case. We can solve them separately:
[8]:
r1 = fanno_solver('friction_sub', [1.06906031, 0])
r2 = fanno_solver('friction_super', [0.3049965, 0.63306493])
print_fanno(*r1)
print_fanno(*r2)
M [0.5 1. ]
P/P* [2.13808993 1. ]
rho/rho* [1.87082869 1. ]
T/T* [1.14285714 1. ]
P0/P0* [1.33984375 1. ]
U/U* [0.53452248 1. ]
4fL*/D [1.06906031 0. ]
(s*-s)/R [0.292553 0. ]
M [1.99999999 3.99999998]
P/P* [0.40824829 0.13363062]
rho/rho* [0.61237244 0.46770717]
T/T* [0.66666667 0.28571429]
P0/P0* [ 1.68749999 10.71874979]
U/U* [1.63299316 2.13808993]
4fL*/D [0.3049965 0.63306493]
(s*-s)/R [0.52324814 2.37199452]
[ ]: