Create Python examples with units

These examples use the units module to implement the existing examples,
but demonstrate how alternative units can be used for the input and
output quantities.

interfaces/cython/cantera/examples/thermo/isentropic_units.py

@@ -0,0 +1,65 @@
+"""
+Isentropic, adiabatic flow example - calculate area ratio vs. Mach number curve
+
+Requires: cantera >= 2.5.0, matplotlib >= 2.0
+"""
+
+import cantera.units as ct
+import math
+import numpy as np
+ct.units.default_format = ".2F~P"
+label_string = "area ratio\tMach number\ttemperature\tpressure ratio"
+output_string = "{0:.2E~P}\t{1}            {2}\t{3:.2E~P}"
+
+
+def soundspeed(gas):
+    """The speed of sound. Assumes an ideal gas."""
+
+    gamma = gas.cp / gas.cv
+    return np.sqrt(gamma * ct.units.molar_gas_constant
+                     * gas.T / gas.mean_molecular_weight).to("m/s")
+
+
+def isentropic(gas=None):
+    """
+    In this example, the area ratio vs. Mach number curve is computed. If a gas
+    object is supplied, it will be used for the calculations, with the
+    stagnation state given by the input gas state. Otherwise, the calculations
+    will be done for a 10:1 hydrogen/nitrogen mixture with stagnation T0 = 1700.33
+    degrees Fahrenheit, P0 = 10 atm.
+    """
+    if gas is None:
+        gas = ct.Solution('gri30.yaml')
+        gas.TPX = 2160 * ct.units.degR, 10.0 * ct.units.atm, 'H2:1,N2:0.1'
+
+    # get the stagnation state parameters
+    s0 = gas.s
+    h0 = gas.h
+    p0 = gas.P
+
+    mdot = 1 * ct.units.kg / ct.units.s  # arbitrary
+    amin = 1.e14 * ct.units.m**2
+
+    data = []
+
+    # compute values for a range of pressure ratios
+    p_range = np.logspace(-3, 0, 200) * p0
+    for p in p_range:
+
+        # set the state using (p,s0)
+        gas.SP = s0, p
+
+        v = np.sqrt(2.0*(h0 - gas.h)).to("m/s")     # h + V^2/2 = h0
+        area = mdot/(gas.density*v)    # rho*v*A = constant
+        amin = min(amin, area)
+        data.append([area, v/soundspeed(gas), gas.T, p/p0])
+
+    return data, amin
+
+
+if __name__ == "__main__":
+    print(__doc__)
+    data, amin = isentropic()
+    print(label_string)
+    for row in data:
+        print(output_string.format(row[0] / amin, row[1], row[2], row[3]))

interfaces/cython/cantera/examples/thermo/rankine_units.py

@@ -0,0 +1,77 @@
+"""
+A Rankine vapor power cycle
+
+Requires: Cantera >= 2.5.0
+"""
+
+import cantera.units as ct
+
+# parameters
+eta_pump = 0.6 * ct.units.dimensionless     # pump isentropic efficiency
+eta_turbine = 0.8 * ct.units.dimensionless # turbine isentropic efficiency
+p_max = 116.03 * ct.units.psi       # maximum pressure
+
+
+def pump(fluid, p_final, eta):
+    """Adiabatically pump a fluid to pressure p_final, using
+    a pump with isentropic efficiency eta."""
+    h0 = fluid.h
+    s0 = fluid.s
+    fluid.SP = s0, p_final
+    h1s = fluid.h
+    isentropic_work = h1s - h0
+    actual_work = isentropic_work / eta
+    h1 = h0 + actual_work
+    fluid.HP = h1, p_final
+    return actual_work
+
+
+def expand(fluid, p_final, eta):
+    """Adiabatically expand a fluid to pressure p_final, using
+    a turbine with isentropic efficiency eta."""
+    h0 = fluid.h
+    s0 = fluid.s
+    fluid.SP =s0, p_final
+    h1s = fluid.h
+    isentropic_work = h0 - h1s
+    actual_work = isentropic_work * eta
+    h1 = h0 - actual_work
+    fluid.HP = h1, p_final
+    return actual_work
+
+
+def printState(n, fluid):
+    print('\n***************** State {0} ******************'.format(n))
+    print(fluid.report())
+
+
+if __name__ == '__main__':
+    # create an object representing water
+    w = ct.Water()
+
+    # start with saturated liquid water at 80.33 degrees Fahrenheit
+    w.TQ = 540 * ct.units.degR, 0.0 * ct.units.dimensionless
+    h1 = w.h
+    p1 = w.P
+    printState(1, w)
+
+    # pump it adiabatically to p_max
+    pump_work = pump(w, p_max, eta_pump)
+    h2 = w.h
+    printState(2, w)
+
+    # heat it at constant pressure until it reaches the saturated vapor state
+    # at this pressure
+    w.PQ = p_max, 1.0 * ct.units.dimensionless
+    h3 = w.h
+    heat_added = h3 - h2
+    printState(3, w)
+
+    # expand back to p1
+    turbine_work = expand(w, p1, eta_turbine)
+    printState(4, w)
+
+    # efficiency
+    eff = (turbine_work - pump_work)/heat_added
+
+    print('efficiency = ', eff)

interfaces/cython/cantera/examples/thermo/sound_speed.py

@@ -6,6 +6,7 @@
 
 import cantera as ct
 import math
+import numpy as np
 
 
 def equilSoundSpeeds(gas, rtol=1.0e-6, max_iter=5000):
@@ -52,7 +53,7 @@ def equilSoundSpeeds(gas, rtol=1.0e-6, max_iter=5000):
 if __name__ == "__main__":
     gas = ct.Solution('gri30.yaml')
     gas.X = 'CH4:1.00, O2:2.0, N2:7.52'
-    for n in range(27):
-        T = 300.0 + 100.0 * n
+    T_range = np.linspace(300, 2900, 27)
+    for T in T_range:
         gas.TP = T, ct.one_atm
         print(T, equilSoundSpeeds(gas))

interfaces/cython/cantera/examples/thermo/sound_speed_units.py

@@ -0,0 +1,60 @@
+"""
+Compute the "equilibrium" and "frozen" sound speeds for a gas
+
+Requires: cantera >= 2.5.0
+"""
+
+import cantera.units as ct
+import numpy as np
+import math
+
+ct.units.default_format = ".2F~P"
+
+def equilSoundSpeeds(gas, rtol=1.0e-6, max_iter=5000):
+    """
+    Returns a tuple containing the equilibrium and frozen sound speeds for a
+    gas with an equilibrium composition.  The gas is first set to an
+    equilibrium state at the temperature and pressure of the gas, since
+    otherwise the equilibrium sound speed is not defined.
+    """
+
+    # set the gas to equilibrium at its current T and P
+    gas.equilibrate('TP', rtol=rtol, max_iter=max_iter)
+
+    # save properties
+    s0 = gas.s
+    p0 = gas.P
+    r0 = gas.density
+
+    # perturb the pressure
+    p1 = p0*1.0001
+
+    # set the gas to a state with the same entropy and composition but
+    # the perturbed pressure
+    gas.SP = s0, p1
+
+    # frozen sound speed
+    afrozen = np.sqrt((p1 - p0)/(gas.density - r0)).to("ft/s")
+
+    # now equilibrate the gas holding S and P constant
+    gas.equilibrate('SP', rtol=rtol, max_iter=max_iter)
+
+    # equilibrium sound speed
+    aequil = np.sqrt((p1 - p0)/(gas.density - r0)).to("ft/s")
+
+    # compute the frozen sound speed using the ideal gas expression as a check
+    gamma = gas.cp/gas.cv
+    afrozen2 = np.sqrt(gamma * ct.units.molar_gas_constant * gas.T /
+                         gas.mean_molecular_weight).to("ft/s")
+
+    return aequil, afrozen, afrozen2
+
+# test program
+if __name__ == "__main__":
+    gas = ct.Solution('gri30.yaml')
+    gas.X = 'CH4:1.00, O2:2.0, N2:7.52'
+    T_range = np.linspace(80.33, 4760.33, 50) * ct.units.degF
+    print("Temperature      Equilibrium Sound Speed     Frozen Sound Speed      Frozen Sound Speed Check")
+    for T in T_range:
+        gas.TP = T, 1.0 * ct.units.atm
+        print(T, *equilSoundSpeeds(gas), sep = "               ")