Fanno equations in aero.Fanno module, usage example

.gitignore

@@ -1,2 +1,4 @@
 # python compiled files
 *.pyc
+# doxygen output folder
+doc

README.md

@@ -7,6 +7,7 @@ Python packages for compressible flow computations
 * Mach dependent functions for isentropic total pressure, temperature and mass flow
 * local Rankine-Hugoniot shock wave equations
 * conical shock waves
+* Fanno equations for momentum losses in a duct
 * misc: degree based trigo functions, Newton iterative solve, ODE integration
 
 #### Installation & usage

aero/Fanno.py

@@ -0,0 +1,48 @@
+"""@package Fanno
+@brief functions to compute Fanno transformation in subsonic/supersonic ducts
+"""
+
+import numpy as np
+
+def maxFparam_Mach(Mach, gamma=1.4):
+    """ computes maximum Fanno parameter for a given Mach number state
+    This maximum value makes the flow reach a sonic state
+    \f$
+      \left(\frac{1 - M^2}{\gamma M^2}\right) + \left(\frac{\gamma + 1}{2\gamma}\right)\ln\left[\frac{M^2}{\left(\frac{2}{\gamma + 1}\right)\left(1 + \frac{\gamma - 1}{2}M^2\right)}\right]
+    \f$
+    """
+    m2 = np.square(Mach)
+    return (1-m2)/(gamma*m2) + (gamma+1)/(2*gamma)*np.log(m2*(gamma+1)/2/(1+(gamma-1)/2*m2))
+
+# --- ratio to critical/choking state ---
+
+def Ps_Pscri(Mach, gamma=1.4):
+    """ computes Ps/Ps_cri ratio from Mach number
+    """
+    return 1./(Mach*np.sqrt(2/(gamma+1)*(1+(gamma-1)/2*np.square(Mach))))
+
+def Ts_Tscri(Mach, gamma=1.4):
+    """ computes Ts/Ts_cri ratio from Mach number
+    """
+    return (gamma+1)/2/(1+(gamma-1)/2*np.square(Mach))
+
+def Rho_Rhocri(Mach, gamma=1.4):
+    """ computes rho/rho_cri ratio from Mach number
+    """
+    return 1./Mach*np.sqrt(2/(gamma+1)*(1+(gamma-1)/2*np.square(Mach)))
+
+def Pi_Picri(Mach, gamma=1.4):
+    """ computes Pi/Pi_cri ratio from Mach number
+    """
+    return 1./Mach*(2/(gamma+1)*(1+(gamma-1)/2*np.square(Mach)))**((gamma+1)/2/(gamma-1))
+
+def V_Vcri(Mach, gamma=1.4):
+    """ computes V/a_cri ratio from Mach number
+    """
+    return Mach / np.sqrt(2/(gamma+1)*(1+(gamma-1)/2*np.square(Mach)))
+
+def NormdS(Mach, gamma=1.4):
+    """ computes delta S / Cp up to critical point
+    """
+    return np.log(Mach**((gamma-1)/gamma) * (2/(gamma+1)*(1+(gamma-1)/2*np.square(Mach)))**(-(gamma+1)/2/gamma))
+

aero/Propulsion.py

@@ -0,0 +1,15 @@
+"""@package Propulsion
+  thrust flow function
+"""
+
+import Isentropic
+import MassFlow
+
+def ThrustFunction(Mach, gamma=1.4):
+    return MassFlow.Sigma_Mach(Mach, gamma)/Isentropic.PiPs_Mach(Mach, gamma)*(1.+gamma*Mach**2)
+
+def Mach_ThrustFunction(thrust, Mach=2., gamma=1.4):
+    def thrust_of_mach(m):
+        return ThrustFunction(m, gamma)
+    return IterativeSolve.secant_solve(thrust_of_mach, thrust, Mach)
+

examples/fanno-diag.py

@@ -0,0 +1,62 @@
+# -*- coding: utf-8 -*-
+"""
+Plot of Fanno curves in several diagrams
+@author: j.gressier
+"""
+
+import numpy                 as np
+import matplotlib.pyplot     as plt
+import aero.Fanno            as fanno
+
+npoints = 100
+gam     = 1.4
+
+Mmin = 0.1
+Mmax = 4.
+
+Mach = np.log10(np.logspace(Mmin, Mmax, npoints+1))
+Fparam = fanno.maxFparam_Mach(Mach, gam)
+Ts     = fanno.Ts_Tscri(Mach, gam)
+Ps     = fanno.Ps_Pscri(Mach, gam)
+Pi     = fanno.Pi_Picri(Mach, gam)
+V      = fanno.V_Vcri(Mach, gam)
+dS     = fanno.NormdS(Mach, gam)
+
+fig=plt.figure(1, figsize=(10,8))
+fig.suptitle('Ratio to critical state, $\gamma = %.1f$'%gam, fontsize=12, y=0.93)
+plt.plot(Mach, Fparam, 'k--')
+plt.plot(Mach, Ts, 'r-')
+plt.plot(Mach, Ps, '-', color='#000088')
+plt.plot(Mach, Pi, '-', color='#0000ff')
+plt.plot(Mach, V, 'k-', color='#009999')
+plt.legend(['Fanno parameter', '$T_s/T^\star$', '$p_s/p^\star$', '$p_i/p_i^\star$', '$V/V^\star$'],
+           loc='upper left',prop={'size':10})  
+plt.axis([Mmin, Mmax, 0., 4.])
+plt.xlabel('Mach number', fontsize=10)
+#plt.ylabel('shock angle $\sigma$', fontsize=10)
+#plt.minorticks_on()
+plt.grid(which='major', linestyle=':', alpha=0.5)
+#plt.grid(which='minor', linestyle=':', alpha=0.5)
+fig.savefig('Fanno-ratio.pdf', bbox_inches='tight')
+#show()
+
+fig=plt.figure(2, figsize=(10,8))
+fig.suptitle('Fanno curve in T/S diagram, $\gamma = %.1f$'%gam, fontsize=12, y=0.93)
+plt.plot(dS, Ts, 'k')
+#==============================================================================
+# Mlab   = np.append(np.linspace(0.1, 1, 10), np.linspace(1, 4, 7))
+# Tslab  = fanno.Ts_Tscri(Mlab, gam)
+# dSlab  = fanno.NormdS(Mlab, gam)
+# plt.plot(dSlab, Tslab, 'ro')
+# for i in range(len(Mlab)):
+#     plt.text(dSlab[i], Tslab[i], '  M=%.2g'%Mlab[i], horizontalalignment='left', verticalalignment='center',
+#              fontsize=8, rotation=-30*np.sign(Mlab[i]-1))
+#==============================================================================
+#plt.axis([Mmin, Mmax, 0., 4.])
+plt.xlabel('$\Delta S/C_p$', fontsize=10)
+plt.ylabel('$h/C_p$', fontsize=10)
+#plt.minorticks_on()
+plt.grid(which='major', linestyle=':', alpha=0.5)
+#plt.grid(which='minor', linestyle=':', alpha=0.5)
+fig.savefig('Fanno-TS.pdf', bbox_inches='tight')
+#show()