.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_shocked_nozzle.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_shocked_nozzle.py: Shocked nozzle design ===================== This example is the inverse of the :doc:`nozzle analysis `: instead of prescribing the area variation and solving for the Mach number, we prescribe the Mach number distribution through a converging--diverging nozzle containing a normal shock and recover the area variation that produces it. The Mach number rises linearly from the inlet, through sonic at the geometric throat, to a supersonic value just upstream of a normal shock. The Rankine--Hugoniot jump conditions return the subsonic state downstream, from which the Mach number falls linearly back to the inlet value. Each isentropic stretch is solved with the same Picard sweep as the analysis example, but driven by the prescribed Mach number (:math:`V = \mathit{Ma}\,a` on the current sound speed) rather than by mass conservation. The shock is the only place an algebraic solve is needed; everything else stays a fixed-point iteration on the equation of state, so the method remains working-fluid independent. .. GENERATED FROM PYTHON SOURCE LINES 20-27 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import ember.block import ember.fluid .. GENERATED FROM PYTHON SOURCE LINES 28-34 Stagnation reservoir -------------------- As before, the upstream reservoir fixes the stagnation enthalpy and the entropy of the (isentropic) flow ahead of the shock, obtained with bare :class:`ember.fluid` calls. .. GENERATED FROM PYTHON SOURCE LINES 34-44 .. code-block:: Python fluid = ember.fluid.PerfectFluid(cp=1005.0, gamma=1.4, mu=1.8e-5, Pr=1.0) Po = 1e5 # Stagnation pressure [Pa] To = 300.0 # Stagnation temperature [K] rhoo, uo = fluid.set_P_T(Po, To) ho = fluid.get_h(rhoo, uo) # Stagnation enthalpy (conserved across the shock) s_up = fluid.get_s(rhoo, uo) # Entropy upstream of the shock .. GENERATED FROM PYTHON SOURCE LINES 45-51 Prescribed Mach number, upstream of the shock --------------------------------------------- The duct is discretised into ``ni`` stations with the shock at mid-length. Upstream, the Mach number rises linearly from the inlet value, passing through unity at the geometric throat, to a supersonic value just before the shock. .. GENERATED FROM PYTHON SOURCE LINES 51-61 .. code-block:: Python ni = 201 x = np.linspace(0.0, 1.0, ni) i_shock = ni // 2 Ma_in = 0.3 # Inlet Mach number Ma_pre = 2.0 # Supersonic Mach number just upstream of the shock Ma_up = np.interp(x[: i_shock + 1], [0.0, x[i_shock]], [Ma_in, Ma_pre]) .. GENERATED FROM PYTHON SOURCE LINES 62-72 Isentropic Picard solver ------------------------ A small helper converges the static state of an isentropic stretch to a prescribed Mach distribution. Each sweep sets the velocity to :math:`\mathit{Ma}` times the *current* sound speed, applies the energy equation :math:`h = h_0 - V^2/2`, and refreshes the state -- the fixed point is exactly the state whose Mach number matches the target, found with no table or per-point inversion. The velocity update is under-relaxed for stability. .. GENERATED FROM PYTHON SOURCE LINES 72-96 .. code-block:: Python def axial_velocity(V): """Pack a scalar axial velocity field into a (n, 3) polar velocity array.""" zero = np.zeros_like(V) return np.stack([V, zero, zero], axis=-1) def solve_isentrope(Ma, s, relax=0.5): """Converge an isentropic Block flow field to a prescribed Mach number.""" n = Ma.size block = ember.block.Block(shape=(n,)).set_fluid(fluid) block.set_h_s(ho * np.ones(n), s).set_Vxrt(axial_velocity(np.zeros(n))) for _ in range(500): V_prev = block.V V = V_prev + relax * (Ma * block.a - V_prev) # Drive velocity from Mach block.set_h_s(ho - 0.5 * V**2, s).set_Vxrt(axial_velocity(V)) if np.max(np.abs(V - V_prev)) < 1e-4: break return block upstream = solve_isentrope(Ma_up, s_up) .. GENERATED FROM PYTHON SOURCE LINES 97-106 Normal shock (Rankine--Hugoniot) -------------------------------- The shock conserves mass, momentum and energy. Taking the upstream state just before the shock, the downstream state is the non-trivial root of the jump conditions. Parametrising by the downstream density :math:`\rho_2`, continuity gives :math:`V_2 = \rho_1 V_1 / \rho_2`, momentum gives :math:`p_2`, and energy gives :math:`h_2`; the equation of state must then reproduce :math:`\rho_2`. A 1-D root find on that residual returns the compressed (subsonic) branch. .. GENERATED FROM PYTHON SOURCE LINES 106-146 .. code-block:: Python rho1, V1, P1 = upstream.rho[-1], upstream.V[-1], upstream.P[-1] mass_flux = rho1 * V1 # rho V, conserved across the shock impulse = P1 + rho1 * V1**2 # p + rho V^2, conserved across the shock def shock_residual(rho2): V2 = mass_flux / rho2 # Continuity P2 = impulse - mass_flux * V2 # Momentum h2 = ho - 0.5 * V2**2 # Energy (ho conserved) rho_eos, _ = fluid.set_P_h(P2, h2) return rho_eos - rho2 def bisect(fun, lo, hi, n=60): """Root of a monotonic residual by bisection (avoids a scipy dependency). A fixed number of halvings is used rather than an absolute tolerance: the block data is single precision, so the interval cannot shrink below the float32 spacing and a tolerance test could otherwise never be satisfied. """ f_lo = fun(lo) for _ in range(n): mid = 0.5 * (lo + hi) if (fun(mid) > 0.0) == (f_lo > 0.0): lo = mid else: hi = mid return 0.5 * (lo + hi) # Bracket above rho1: a shock compresses the gas (rho2 > rho1). rho2 = bisect(shock_residual, rho1 * (1.0 + 1e-4), rho1 * 20.0) V2 = mass_flux / rho2 P2 = impulse - mass_flux * V2 h2 = ho - 0.5 * V2**2 _, u2 = fluid.set_P_h(P2, h2) s_down = fluid.get_s(rho2, u2) # Entropy rises across the shock Ma_down = V2 / fluid.get_a(rho2, u2) # Subsonic post-shock Mach number .. GENERATED FROM PYTHON SOURCE LINES 147-154 Prescribed Mach number, downstream of the shock ----------------------------------------------- Downstream the flow is isentropic again, but on the higher entropy ``s_down`` (the stagnation pressure has dropped while the stagnation enthalpy is unchanged). The Mach number falls linearly from the post-shock value back to the inlet value, and the same solver gives the state. .. GENERATED FROM PYTHON SOURCE LINES 154-161 .. code-block:: Python Ma_dn = np.interp(x[i_shock:], [x[i_shock], 1.0], [Ma_down, Ma_in]) downstream = solve_isentrope(Ma_dn, s_down) print(f"Shock: Ma {Ma_pre:.2f} -> {Ma_down:.3f}") print(f"Stagnation pressure ratio across shock: {downstream.Po[0] / Po:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Shock: Ma 2.00 -> 0.577 Stagnation pressure ratio across shock: 0.721 .. GENERATED FROM PYTHON SOURCE LINES 162-172 Area variation from mass conservation ------------------------------------- With the states known everywhere, the area follows algebraically from continuity :math:`\rho V A = \dot{m}`: normalising by the inlet, the area ratio is the inlet mass flux divided by the local mass flux. The mass flux is continuous across the shock, so the area is too -- even though the Mach number jumps. Note that the exit area exceeds the inlet area despite ``Ma`` returning to its inlet value, because the entropy rise leaves the two stations on different isentropes. .. GENERATED FROM PYTHON SOURCE LINES 172-181 .. code-block:: Python G_in = (upstream.rho * upstream.V)[0] A_up = G_in / (upstream.rho * upstream.V) A_dn = G_in / (downstream.rho * downstream.V) i_throat = np.argmin(A_up) # Geometric throat = minimum area (sonic point) print(f"Throat area ratio A/A_in = {A_up[i_throat]:.3f}") print(f"Exit area ratio A/A_in = {A_dn[-1]:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Throat area ratio A/A_in = 0.491 Exit area ratio A/A_in = 1.387 .. GENERATED FROM PYTHON SOURCE LINES 182-188 Result ------ The prescribed Mach distribution (with its shock discontinuity) and the recovered area variation, drawn against axial position. The throat and shock locations are marked. .. GENERATED FROM PYTHON SOURCE LINES 188-206 .. code-block:: Python fig, axs = plt.subplots(2, 1, sharex=True, figsize=(6, 5.5)) axs[0].plot(x[: i_shock + 1], upstream.Ma, "-") axs[0].plot(x[i_shock:], downstream.Ma, "-", color="C0") axs[0].set_ylabel("Mach number, $\\mathit{Ma}$ [-]") axs[1].plot(x[: i_shock + 1], A_up, "-") axs[1].plot(x[i_shock:], A_dn, "-", color="C0") axs[1].set_ylabel("Area ratio, $A/A_\\mathrm{in}$ [-]") axs[1].set_xlabel("Axial position, $x$ [-]") for ax in axs: ax.axvline(x[i_throat], color="0.7", ls="--", lw=1.0) # Throat ax.axvline(x[i_shock], color="C3", ls=":", lw=1.0) # Shock fig.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_shocked_nozzle_001.png :alt: plot shocked nozzle :srcset: /auto_examples/images/sphx_glr_plot_shocked_nozzle_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.158 seconds) .. _sphx_glr_download_auto_examples_plot_shocked_nozzle.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_shocked_nozzle.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_shocked_nozzle.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_shocked_nozzle.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_