Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Flows with heat transfer (Rayleigh flows)

Notebook Cell
%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np
from scipy.optimize import root_scalar

# Pint gives us some helpful unit conversion
from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity # We will use this to construct quantities (value + unit)
Notebook Cell
# these lines are only for helping improve the display
import matplotlib_inline.backend_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('pdf', 'png')
plt.rcParams['figure.dpi']= 300
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['mathtext.fontset'] = 'cm'

Heat transfer is the third main factor that can affect a compressible flow, after area change and friction.

Rayleigh flow is the specific case of frictionless flow in a constant-area duct, with heat transfer. Rayleigh flow applies to constant-area heat exchangers and combustion chambers where the entropy changes due to heat transfer are significantly larger than those due to friction:

dsedsi  ,d s_e \gg d s_i \;,

so we can neglect frictional effects and say dsdse ds \approx ds_e .

Rayleigh flow theory

Control volume for Rayleigh flow

Figure 8:Differential control volume for frictionless flow in a constant-area duct with heat transfer.

Figure 8 shows a control volume in a constant-area duct, where the flow starts out with velocity VV, density ρ\rho, and pressure pp. A small amount of heat δq\delta q is added to the flow. The width of the differential control volume is dxdx and the cross-sectional area of the duct is AA.

For a general fluid, we can apply conservation of mass:

ρV=(ρ+dρ)(V+dV)ρV=G=constant  .\begin{gather*} \rho V = (\rho + d\rho) (V + dV) \rightarrow \rho V = G = \text{constant} \;. \end{gather*}

Applying conservation of momentum gives

Fx=m˙(Vout,xVin,x)pA(p+dp)A=(ρAV)[(V+dV)V]Adp=ρAVdVdp+ρVdV=0  ,\begin{gather*} \sum F_x = \dot{m} \left( V_{\text{out}, x} - V_{\text{in}, x} \right) \\ pA - (p + dp) A = (\rho A V) \left[ (V + dV) - V \right] \\ - A dp = \rho A V dV \\ dp + \rho V dV = 0 \;, \end{gather*}

which we can integrate through the flow (since ρV=\rho V = constant) to get

p+ρV2=constantp + \rho V^2 = \text{constant}

Next, applying conservation of energy:

δqδws=[(h+dh)+(V+dV)22](h+V22)δq=dh+VdV  .\begin{gather*} \delta q - \delta w_s = \left[ (h+dh) + \frac{(V+dV)^2}{2} \right] - \left( h + \frac{V^2}{2} \right) \\ \delta q = dh + V dV \;. \end{gather*}

For an ideal gas, where dh=cpdT dh = c_p dT , we can express the above as

δq=cpdT+VdV  .\delta q = c_p dT + V dV \;.

Recall the definitions of stagnation enthalpy and temperature:

ht=cpTt=h+V22=cpT+V22  .h_t = c_p T_t = h + \frac{V^2}{2} = c_p T + \frac{V^2}{2} \;.

If we differentiate, we get

dht=cpdTt=cpT+VdVd h_t = c_p d T_t = c_p T + V dV

and see that

δq=dht=cpdTt  .\delta q = d h_t = c_p d T_t \;.

We can see that stagnation enthlapy is not constant.

We can use these relationships to develop a TTss or hhss diagram of a Rayleigh line, based on the property relationship

Tds=dhvdpds=cpdTTRdpp  ,\begin{gather*} T ds = dh - v dp \\ ds = c_p \frac{dT}{T} - R \frac{dp}{p} \;, \end{gather*}

which we can integrate from some reference state 1 to some general state:

ss1cp=logTT1γ1γlogpp1  .\frac{s-s_1}{c_p} = \log \frac{T}{T_1} - \frac{\gamma-1}{\gamma} \log \frac{p}{p_1} \;.

To express the entropy difference (ss1)/cp (s - s_1)/c_p only in terms of relative temperature T/T1 T / T_1 and the reference Mach number M1 M_1 we need to do a bit of work, but we can obtain

ss1cp=logTT1γ1γlog[(1+γM12)±(1+γM12)24γM12(T/T1)2]  ,\frac{s - s_1}{c_p} = \log \frac{T}{T_1} - \frac{\gamma - 1}{\gamma} \log \left[ \frac{(1 + \gamma M_1^2) \pm \sqrt{ \left(1 + \gamma M_1^2 \right)^2 - 4 \gamma M_1^2 (T/T_1)}}{2} \right] \;,

where the ++ sign in the equation is used for the upper branch of the Rayleigh line to the point of maximum temperature, and the - sign in the equation is used for the lower part of the curve. Given a reference Mach number M1M_1 and specific heat ratio γ\gamma, we can use Equation (215) to compute a Rayleigh line on a TTss diagram, using relative temperature T/T1T / T_1 and entropy (ss1)/cp (s - s_1)/c_p .

gamma = 1.3
M1 = 3.0
Tt_T1 = 1 + 0.5*(gamma-1)*M1**2

max_temp_ratio = (1+gamma*M1**2)**2 / (4*gamma*M1**2)

T_T1 = np.linspace(1.0, max_temp_ratio, 300)
delta_s_cp_upper = (
    np.log(T_T1) - (gamma-1)*np.log(
        (1+gamma*M1**2 + np.sqrt((1+gamma*M1**2)**2 - 4*gamma*T_T1*M1**2))/2
        ) / gamma
    )
delta_s_cp_lower = (
    np.log(T_T1) - (gamma-1)*np.log(
        (1+gamma*M1**2 - np.sqrt((1+gamma*M1**2)**2 - 4*gamma*T_T1*M1**2))/2
        ) / gamma
    )

fig, ax = plt.subplots()

ax.plot(
    np.hstack((delta_s_cp_lower, delta_s_cp_upper[::-1])),
    np.hstack((T_T1, T_T1[::-1]))
    )
ax.set_title(r'Rayleigh line ($\gamma = 1.3$)')
ax.set_xlabel(r'$(s - s_1) / c_p$')
ax.set_ylabel(r'$T / T_1$')
ax.set_xlim([0, 0.9])
ax.grid(True)
#plt.savefig('rayleigh-line.png')
plt.tight_layout()
plt.show()
#plt.close()
Loading...

As with Fanno flows, we see that there is a limiting point of maximum entropy. However, unlike with Fanno flows, we can move both to the left and right along a Rayleigh line, because heat can be added or removed. Adding heat means positive entropy and moving towards the limiting point, and removing heat means entropy decreases and moving away from the limiting point.

Identifying the state associated with this limiting point follows a similar procedure as with Fanno flows. From conservation of mass, we have G=ρV= G = \rho V = constant. From conservation of momentum given by Equation (207), we have

p+ρV2=G2ρ=constantdp+G2(dρρ2)=0dpdρ=G2ρ2=V2\begin{gather*} p + \rho V^2 = \frac{G^2}{\rho} = \text{constant} \\ dp + G^2 \left(- \frac{d\rho}{\rho^2} \right) = 0 \\ \frac{dp}{d\rho} = \frac{G^2}{\rho^2} = V^2 \\ \end{gather*}

This relationship applies to any fluid, between two differentially separated states on a Rayleigh line. If we apply this at two points right around the limiting point of maximum entropy, then entropy is constant and ds=0 ds = 0 :

V2=(dpdρ)limiting point=(pρ)sV2=a2\begin{gather*} V^2 = \left( \frac{dp}{d\rho} \right)_{\text{limiting point}} = \left( \frac{\partial p}{\partial \rho} \right)_{s} \\ \therefore V^2 = a^2 \end{gather*}

and the velocity is sonic at the limiting point.

Working relations for Rayleigh flows

Based on the analysis above, we can establish working relations for properties in a Rayleigh flow as a function of Mach number.

Static properties

From conservation of mometum, we already have an expression for pressure given by Equation (217), which we can express between two points as

p2p1=1+γM121+γM22  .\frac{p_2}{p_1} = \frac{1 + \gamma M_1^2}{1 + \gamma M_2^2} \;.

We also already developed a relationship between temperature at two points in the flow in Equation (222), which we can express in general as

T2T1=M22M12(1+γM121+γM22)2\frac{T_2}{T_1} = \frac{M_2^2}{M_1^2} \left( \frac{1 + \gamma M_1^2}{1 + \gamma M_2^2} \right)^2

To get a relationship for density, we can use conservation of mass and the speed of sound for an ideal gas:

ρ1V1=ρ2V2ρ2ρ1=V1V2=M1γRT1M2γRT2=M1M2T1T2\begin{gather*} \rho_1 V_1 = \rho_2 V_2 \\ \frac{\rho_2}{\rho_1} = \frac{V_1}{V_2} = \frac{M_1 \sqrt{\gamma R T_1}}{M_2 \sqrt{\gamma R T_2}} = \frac{M_1}{M_2} \sqrt{\frac{T_1}{T_2}} \end{gather*}

and then substitute Equation (229) to get

ρ2ρ1=M12M22(1+γM221+γM22)  .\frac{\rho_2}{\rho_1} = \frac{M_1^2}{M_2^2} \left( \frac{ 1 + \gamma M_2^2}{1 + \gamma M_2^2 } \right) \;.

Stagnation properties

We know that the stagnation temperature does not remain constant in a Rayleigh flow, so we need an expression for stagnation temperature ratio. Using the definition of stagnation temperature

Tt=T(1+γ12M2)  ,T_t = T \left( 1 + \frac{\gamma - 1}{2} M^2 \right) \;,

we can write this for locations 1 and 2 and take the ratio to get

Tt2Tt1=T2T1(1+γ12M221+γ12M12)  ,\frac{T_{t2}}{T_{t1}} = \frac{T_2}{T_1} \left( \frac{1 + \frac{\gamma-1}{2} M_2^2}{1 + \frac{\gamma-1}{2} M_1^2} \right) \;,

which we can combine with the static temperature ratio given by Equation (229) to get

Tt2Tt1=M22M12(1+γM121+γM22)2(1+γ12M221+γ12M12)  .\frac{T_{t2}}{T_{t1}} = \frac{M_2^2}{M_1^2} \left( \frac{1 + \gamma M_1^2}{1 + \gamma M_2^2} \right)^2 \left( \frac{1 + \frac{\gamma-1}{2} M_2^2}{1 + \frac{\gamma-1}{2} M_1^2} \right) \;.

Similarly, we can get the stagnation pressure ratio by taking the definition of stagnation pressure

pt=p(1+γ12M2)γ/(γ1)p_t = p \left( 1 + \frac{\gamma-1}{2} M^2 \right)^{\gamma/(\gamma-1)}

and taking the ratio between two locations:

pt2pt1=p2p1(1+γ12M221+γ12M12)γ/(γ1)  .\frac{p_{t2}}{p_{t1}} = \frac{p_2}{p_1} \left( \frac{1 + \frac{\gamma-1}{2} M_2^2}{1 + \frac{\gamma-1}{2} M_1^2} \right)^{\gamma / (\gamma - 1)} \;.

Then, substitute the static pressure ratio in Equation (228):

pt2pt1=1+γM121+γM22(1+γ12M221+γ12M12)γ/(γ1)  .\frac{p_{t2}}{p_{t1}} = \frac{1 + \gamma M_1^2}{1 + \gamma M_2^2} \left( \frac{1 + \frac{\gamma-1}{2} M_2^2}{1 + \frac{\gamma-1}{2} M_1^2} \right)^{\gamma / (\gamma - 1)} \;.

We can write functions to evaluate these property ratios:

def get_rayleigh_pressure_ratio(mach_1, mach_2, gamma=1.4):
    '''Return p2/p1 for Rayleigh flow'''
    return ((1 + gamma*mach_1**2) / (1 + gamma*mach_2**2))

def get_rayleigh_temperature_ratio(mach_1, mach_2, gamma=1.4):
    '''Return T2/T1 for Rayleigh flow'''
    return (
        ((1 + gamma*mach_1**2) / (1 + gamma*mach_2**2))**2 *
        (mach_2**2 / mach_1**2)
        )
            
def get_rayleigh_density_ratio(mach_1, mach_2, gamma=1.4):
    '''Return rho2/rho1 for Rayleigh flow'''
    return (
        ((1 + gamma*mach_2**2) / (1 + gamma*mach_1**2)) *
        (mach_1**2 / mach_2**2)
        )
            
def get_rayleigh_stag_temperature_ratio(mach_1, mach_2, gamma=1.4):
    '''Return Tt2/Tt1 for Rayleigh flow'''
    return (
        ((1 + gamma*mach_1**2) / (1 + gamma*mach_2**2))**2 * 
            (mach_2 / mach_1)**2 * (
            (1 + 0.5*(gamma-1)*mach_2**2) /
            (1 + 0.5*(gamma-1)*mach_1**2)
            )
        )

def get_rayleigh_stag_pressure_ratio(mach_1, mach_2, gamma=1.4):
    '''Return pt2/pt1 for Rayleigh flow'''
    return (
        ((1 + gamma*mach_1**2) / (1 + gamma*mach_2**2)) * (
            (1 + 0.5*(gamma-1)*mach_2**2) /
            (1 + 0.5*(gamma-1)*mach_1**2)
            )**(gamma / (gamma - 1))
        )

Energy

The property relationships given by Equations (228), (229), (231), (234), and (237) all require knowledge of the downstream Mach number. If we do not know those conditions, we need to find them by considering conservation of energy and the heat transfer in the flow.

If we apply conservation of energy between two locations in the flow, with some heat transfer qq, recalling that h=cpT h = c_p T and ht=cpTt h_t = c_p T_t , we get

ht1+q=ht2cpTt1+q=cpTt2\begin{gather*} h_{t1} + q = h_{t2} \\ c_p T_{t1} + q = c_p T_{t2} \end{gather*}

and thus

q=cp(Tt2Tt1)  .q = c_p \left( T_{t2} - T_{t1} \right) \;.

Given an amount of heat transfer, this relationship can be used to relate the two locations in the flow.

Rayleigh sonic reference state

The * reference state for a Rayleigh flow is the state that would exist if the flow continued until the Mach number is 1.0, through additional heat transfer. Since this is the limiting point on a Rayleigh line, all states in a given Rayleigh flow share the same * reference state. This can be a useful concept to solve Rayleigh flow problems, because we can easily compute property ratios between the flow at a given location and the reference point.

If we take our property relationships and apply them between an arbitrary point in the flow system and the Rayleigh * reference state, we can get

pp=f(M,γ)=1+γ1+γM2\frac{p}{p^*} = f(M, \gamma) = \frac{1 + \gamma}{1 + \gamma M^2}
TT=f(M,γ)=M2(1+γ)2(1+γM2)2\frac{T}{T^*} = f(M, \gamma) = \frac{M^2 (1 + \gamma)^2}{(1 + \gamma M^2)^2}
ρρ=f(M,γ)=1+γM2(1+γ)M2\frac{\rho}{\rho^*} = f(M, \gamma) = \frac{1 + \gamma M^2}{(1 + \gamma) M^2}
TtTt=f(M,γ)=2(1+γ)M2(1+γM2)2(1+γ12M2)\frac{T_t}{T_t^*} = f(M, \gamma) = \frac{2(1 + \gamma) M^2}{(1 + \gamma M^2)^2} \left( 1 + \frac{\gamma-1}{2} M^2 \right)
ptpt=f(M,γ)=1+γ1+γM2(1+γ12M2(γ+1)/2)γ/(γ1)\frac{p_t}{p_t^*} = f(M, \gamma) = \frac{1 + \gamma}{1 + \gamma M^2} \left( \frac{1 + \frac{\gamma-1}{2} M^2}{(\gamma + 1)/2} \right)^{\gamma/(\gamma-1)}
def get_rayleigh_reference_pressure(mach, gamma=1.4):
    '''Return p/p* for Rayleigh flow'''
    return ((1 + gamma) / (1 + gamma*mach**2))

def get_rayleigh_reference_temperature(mach, gamma=1.4):
    '''Return T/T* for Rayleigh flow'''
    return (
        mach**2 * (1 + gamma)**2 / 
        (1 + gamma*mach**2)**2
        )
            
def get_rayleigh_reference_density(mach, gamma=1.4):
    '''Return rho/rho* for Rayleigh flow'''
    return ((1 + gamma*mach**2) / ((1 + gamma) * mach**2))
            
def get_rayleigh_reference_stag_temperature(mach, gamma=1.4):
    '''Return Tt/Tt* for Rayleigh flow'''
    return (
        2*(1 + gamma)*mach**2 * 
        (1 + 0.5*(gamma - 1)*mach**2) / (1 + gamma*mach**2)**2
        )

def get_rayleigh_reference_stag_pressure(mach, gamma=1.4):
    '''Return pt/pt* for Rayleigh flow'''
    return (
        (1 + gamma) * (
            (1 + 0.5*(gamma-1)*mach**2) / (0.5*(gamma+1))
            )**(gamma / (gamma - 1)) /
            (1 + gamma*mach**2)
        )

Example: find heat transfer direction

Consider a constant-area duct with supersonic flow of air, where at one location M1=1.5 M_1 = 1.5 and p1=10 p_1 = 10 bar, and at a downstream location M2=3.0 M_2 = 3.0 . Find the pressure at location 2 and the direction of heat transfer.

We can use two different strategies to find the downstream pressure p2p_2. Using the reference state, which is shared between the two locations, we can write

p2=p2ppp1p1p_2 = \frac{p_2}{p^*} \frac{p^*}{p_1} p_1

and use the Mach numbers at each location to find the two ratios.

gamma = 1.4
M1 = 1.5
M2 = 3.0

p1 = Q_(10, 'bar')

p1_pstar = get_rayleigh_reference_pressure(M1, gamma)
p2_pstar = get_rayleigh_reference_pressure(M2, gamma)
p2 = p2_pstar * (1 / p1_pstar) * p1
print(f'p2 = {p2: .2f}')
p2 =  3.05 bar

Alternatively, we can use the general pressure ratio p2/p1p_2 / p_1 to directly get p2p_2:

p2_p1 = get_rayleigh_pressure_ratio(M1, M2, gamma)
p2 = p2_p1 * p1
print(f'p2 = {p2: .2f}')
p2 =  3.05 bar

To find the direction of heat transfer, we can calculate the stagnation temperature ratio Tt2Tt1 \frac{T_{t2}}{T_{t1}} . We can also do this using either approach.

Using the reference state,

Tt2Tt1=Tt2TtTtTt1\frac{T_{t2}}{T_{t1}} = \frac{T_{t2}}{T_t^*} \frac{T_t^*}{T_{t1}}
Tt1_Ttstar = get_rayleigh_reference_stag_temperature(M1, gamma)
Tt2_Ttstar = get_rayleigh_reference_stag_temperature(M2, gamma)

Tt2_Tt1 = Tt2_Ttstar / Tt1_Ttstar
print(f'T_t2 / T_t1 = {Tt2_Tt1: .3f}')
T_t2 / T_t1 =  0.719

Or, we can directly calculate it:

Tt2_Tt1 = get_rayleigh_stag_temperature_ratio(M1, M2, gamma)
print(f'T_t2 / T_t1 = {Tt2_Tt1: .3f}')
T_t2 / T_t1 =  0.719

Since the stagnation pressure is decreasing, that means that the flow is cooling.

Example: flow through combustion chamber

Consider flow of air into a combustion chamber, where we represent the effects of combustion via heat addition at 120 kJ/kg. The conditions of the flow entering the chamber are 220 K, 70 kPa, and it has a velocity of 122 m/s. What are the exit conditions?

First, we need to find the Mach number and stagnation temperature at the inlet. We can also identify the change in stagnation temperature using q=cpΔTt q = c_p \Delta T_t , since we know that for air cpc_p = 1000 J/(kg K).

def get_stagnation_temperature_ratio(mach, gamma=1.4):
    '''Calculates T/Tt stagnation relationship.'''
    return (1.0 / (1 + 0.5*(gamma - 1) * mach**2))
gamma = 1.4
R = Q_(287, 'J/(kg*K)')
cp = Q_(1000, 'J/(kg*K)')
q = Q_(120, 'kJ/kg')

T1 = Q_(220, 'K')
p1 = Q_(70, 'kPa')
V1 = Q_(122, 'm/s')

a1 = np.sqrt(gamma * R * T1)
M1 = V1 / a1
print(f'Mach @ 1 = {M1.to_base_units(): .3f~P}')
T1_Tt1 = get_stagnation_temperature_ratio(M1, gamma)
Tt1 = T1 / T1_Tt1

delta_Tt = q / cp
Tt2 = Tt1 + delta_Tt
Mach @ 1 =  0.410

To find the downstream conditions, we have two options for solution approaches.

Solution via reference state

The two locations share the same Rayleigh sonic reference state, so we can write

Tt2Tt=Tt2Tt1Tt1Tt=f(M2,γ)p2=p2ppp1T2=T2TTT1\begin{gather*} \frac{T_{t2}}{T_t^*} = \frac{T_{t2}}{T_{t1}} \frac{T_{t1}}{T_t^*} = f(M_2, \gamma) \\ p_2 = \frac{p_2}{p^*} \frac{p^*}{p_1} \\ T_2 = \frac{T_2}{T^*} \frac{T^*}{T_1} \end{gather*}

To solve for M2 M_2 given Tt2/T T_{t2} / T^* , we need to write a function to use root_scalar():

def solve_rayleigh_reference_stagnation_temperature(mach, Tt_Ttstar, gamma=1.4):
    '''Used to find unknown Mach number given Tt/Tt* and gamma'''
    return (
        Tt_Ttstar - get_rayleigh_reference_stag_temperature(mach, gamma)
        )
Tt1_Tstar = get_rayleigh_reference_stag_temperature(M1, gamma)
T1_Tstar = get_rayleigh_reference_temperature(M1, gamma)
p1_pstar = get_rayleigh_reference_pressure(M1, gamma)

Tt2_Tstar = (Tt2 / Tt1) * Tt1_Tstar

# we know the flow is subsonic
root = root_scalar(
    solve_rayleigh_reference_stagnation_temperature, x0=0.5, x1=0.6,
    args=(Tt2_Tstar, gamma)
    )
M2 = root.root
print(f'Mach @ 2 = {M2.to_base_units(): .3f~P}')

T2_Tstar = get_rayleigh_reference_temperature(M2, gamma)
p2_pstar = get_rayleigh_reference_pressure(M2, gamma)

T2 = T2_Tstar * (1/T1_Tstar) * T1
print(f'T2 = {T2.to("K"): .2f~P}')
p2 = p2_pstar * (1/p1_pstar) * p1
print(f'p2 = {p2: .2f~P}')
Mach @ 2 =  0.616
T2 =  322.90 K
p2 =  56.49 kPa

Solution via direct relationships

Alternatively, we can solve this problem by directly solving the property relationships between two states in a Rayleigh flow. Since we know Tt2T_{t2} and Tt1T_{t1}, we can solve Equation (234) for M2M_2. Once we know both Mach numbers, we can use the property relationships for pressure and temperature.

def solve_rayleigh_stagnation_temperature(M2, M1, Tt2_Tt1, gamma):
    '''Used to find unknown M2 given M1, stagnation temperature ratio, and gamma.'''
    return (
        Tt2_Tt1 - get_rayleigh_stag_temperature_ratio(M1, M2, gamma)
        )
Tt2_Tt1 = Tt2 / Tt1
root = root_scalar(
    solve_rayleigh_stagnation_temperature, 
    x0=0.5, x1=0.7, args=(M1, Tt2_Tt1, gamma)
    )
M2 = root.root
print(f'Mach @ 2 = {M2.to_base_units(): .3f~P}')

p2 = p1 * get_rayleigh_pressure_ratio(M1, M2, gamma)
T2 = T1 * get_rayleigh_temperature_ratio(M1, M2, gamma)
print(f'T2 = {T2.to("K"): .2f~P}')
print(f'p2 = {p2: .2f~P}')
Mach @ 2 =  0.616
T2 =  322.90 K
p2 =  56.49 kPa

Both solution approaches lead to the same answers, as expected.