.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/aerodynamics/compressible/plot_09_fanno_flow.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_aerodynamics_compressible_plot_09_fanno_flow.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_aerodynamics_compressible_plot_09_fanno_flow.py:


Fanno flow: friction chokes a duct (1904)
=========================================

Gino Fanno's 1904 diploma thesis at the ETH Zürich, written under Aurel
Stodola, followed an adiabatic gas along a constant-area pipe with wall
friction. Friction drives the Mach number towards one, from below in
subsonic flow and from above in supersonic flow, while the stagnation
temperature stays fixed and the stagnation pressure falls. The states lie on
a *Fanno line* in the temperature-entropy plane, whose point of maximum
entropy is the sonic state. :func:`~engineeringkit.aerodynamics.fanno_flow`
gives every ratio to that sonic state, and
:func:`~engineeringkit.aerodynamics.fanno_choking_length` the pipe length
that chokes it. This example marches the friction ODE down a pipe and
checks the closed form.

.. GENERATED FROM PYTHON SOURCE LINES 19-27

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import solve_ivp

    from engineeringkit.aerodynamics import fanno_choking_length, fanno_flow

    gamma, D, f = 1.4, 0.05, 0.02








.. GENERATED FROM PYTHON SOURCE LINES 28-34

March the Mach number down the pipe
-----------------------------------

Momentum and energy give
:math:`\frac{dM}{dx} = \frac{f}{2D}\,\frac{\gamma M^3(1 + \tfrac{\gamma-1}{2}M^2)}{1 - M^2}`
with the Darcy friction factor :math:`f`.

.. GENERATED FROM PYTHON SOURCE LINES 34-52

.. code-block:: Python



    def dM_dx(x, M):
        return f / (2 * D) * gamma * M**3 * (1 + (gamma - 1) / 2 * M**2) / (1 - M**2)


    fig, axes = plt.subplots(1, 2, figsize=(10, 4))
    for M_in in (0.3, 0.5, 2.0, 3.0):
        L_star = float(fanno_choking_length(M_in, D, f, gamma))
        sonic = solve_ivp(lambda x, M: dM_dx(x, M), (0, 0.999 * L_star), [M_in], rtol=1e-10, atol=1e-12, dense_output=True)
        x = np.linspace(0, 0.999 * L_star, 200)
        axes[0].plot(x / D, sonic.sol(x)[0], label=f"$M_1$ = {M_in}, $L^*$ = {L_star:.2f} m")
        print(f"M1 = {M_in}: L* = {L_star:.3f} m, Mach after 0.999 L* = {sonic.y[0, -1]:.3f}")
        assert abs(sonic.y[0, -1] - 1) < 0.1
    axes[0].axhline(1, color="gray", lw=0.8)
    axes[0].set(xlabel="distance along pipe $x/D$", ylabel="Mach number", title="Friction drives M towards 1")
    axes[0].legend(fontsize=8)




.. image-sg:: /api/gallery/aerodynamics/compressible/images/sphx_glr_plot_09_fanno_flow_001.png
   :alt: Friction drives M towards 1
   :srcset: /api/gallery/aerodynamics/compressible/images/sphx_glr_plot_09_fanno_flow_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    M1 = 0.3: L* = 52.993 m, Mach after 0.999 L* = 1.000
    M1 = 0.5: L* = 10.691 m, Mach after 0.999 L* = 1.000
    M1 = 2.0: L* = 3.050 m, Mach after 0.999 L* = 1.000
    M1 = 3.0: L* = 5.222 m, Mach after 0.999 L* = 1.000

    <matplotlib.legend.Legend object at 0x7fd5b65bb830>



.. GENERATED FROM PYTHON SOURCE LINES 53-59

The Fanno line
--------------

Entropy relative to the sonic state is
:math:`(s - s^*)/c_p = \ln(T/T^*) - \tfrac{\gamma-1}{\gamma}\ln(p/p^*)`. Its
maximum, where friction can do no more, is at :math:`M = 1`.

.. GENERATED FROM PYTHON SOURCE LINES 59-74

.. code-block:: Python


    M = np.linspace(0.15, 3.5, 400)
    state = fanno_flow(M, gamma)
    entropy = np.log(state.temperature_ratio) - (gamma - 1) / gamma * np.log(state.pressure_ratio)
    axes[1].plot(entropy, state.temperature_ratio)
    axes[1].plot(0, 1, "ko")
    axes[1].set(xlabel="$(s - s^*)/c_p$", ylabel="$T/T^*$", title="Fanno line")
    axes[1].annotate("subsonic", (entropy[20], state.temperature_ratio[20]))
    axes[1].annotate("supersonic", (entropy[-60], state.temperature_ratio[-60]))
    fig.tight_layout()

    assert np.isclose(M[np.argmax(entropy)], 1.0, atol=0.01)
    assert np.isclose(float(fanno_flow(2.0).friction_parameter) / 4, 0.30500, atol=1e-5)

    plt.show()








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.192 seconds)


.. _sphx_glr_download_api_gallery_aerodynamics_compressible_plot_09_fanno_flow.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../../lite/lab/index.html?path=api/gallery/aerodynamics/compressible/plot_09_fanno_flow.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_09_fanno_flow.ipynb <plot_09_fanno_flow.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_09_fanno_flow.py <plot_09_fanno_flow.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_09_fanno_flow.zip <plot_09_fanno_flow.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_