Hugoniot equation for a calorically perfect gas

The Hugoniot equation relates thermodynamic quantities (pressure \(p\), specific volume \(v\) and internal energy \(e\)) across a shock wave.

The general form (shown below) holds for a perfect gas, chemically reacting gas, real gas, etc.

from sympy import *
from sympy_equation import (
    Eqn,
    table_of_expressions,
    process_arguments_of_add,
    divide_term_by_term,
)

e1, e2, p1, p2, v1, v2 = symbols("e1, e2, p1, p2, v1, v2", real=True, positive=True)
h = Eqn(e2 - e1, (p1 + p2) / 2 * (v1 - v2))
h

$\displaystyle - e_{1} + e_{2} = \left(\frac{p_{1}}{2} + \frac{p_{2}}{2}\right) \left(v_{1} - v_{2}\right)$

For a calorically perfect gas, \(e = C_{v} T\) and \(T = p v / R\), where:

• \(C_{v}\): specific heat at constant volume.

• \(T\): temperature.

• \(R\): specific gas constant.

Also note that \(\gamma\) is the ratio of specific heats.

Find a new form of the Hugoniot equation for a calorically perfect gas, which shows the relationship between the pressure ratio and the specific volume ratio.

Let’s define a few new equations:

Cv, R, T1, T2, gamma = symbols("C_v, R, T1, T2, gamma", real=True, positive=True)
eq1 = Eqn(Cv, R / (gamma - 1))
eq2 = Eqn(e1, Cv * T1)
eq3 = Eqn(e2, Cv * T2)
eq4 = Eqn(T1, p1 * v1 / R)
eq5 = Eqn(T2, p2 * v2 / R)
display(eq1, eq2, eq3, eq4, eq5)

$\displaystyle C_{v} = \frac{R}{\gamma - 1}$

$\displaystyle e_{1} = C_{v} T_{1}$

$\displaystyle e_{2} = C_{v} T_{2}$

$\displaystyle T_{1} = \frac{p_{1} v_{1}}{R}$

$\displaystyle T_{2} = \frac{p_{2} v_{2}}{R}$

Substitutes them into the Hugoniot equation:

# multi step substitutions in order to achieve the desired form
e = h.subs(eq2, eq3).subs(eq4, eq5).subs(eq1)
e

$\displaystyle - \frac{p_{1} v_{1}}{\gamma - 1} + \frac{p_{2} v_{2}}{\gamma - 1} = \left(\frac{p_{1}}{2} + \frac{p_{2}}{2}\right) \left(v_{1} - v_{2}\right)$

Let’s move \(\gamma - 1\) from the denominator of the LHS to the numerator of the RHS:

e = e.dolhs.collect(1 / (gamma - 1)) * (gamma - 1)
e

$\displaystyle - p_{1} v_{1} + p_{2} v_{2} = \left(\gamma - 1\right) \left(\frac{p_{1}}{2} + \frac{p_{2}}{2}\right) \left(v_{1} - v_{2}\right)$

Now we can start working with pressure ratios and specific volume ratios. First, let’s expand the equation so it is easier to achieve the next steps:

e = e.expand()
e

$\displaystyle - p_{1} v_{1} + p_{2} v_{2} = \frac{\gamma p_{1} v_{1}}{2} - \frac{\gamma p_{1} v_{2}}{2} + \frac{\gamma p_{2} v_{1}}{2} - \frac{\gamma p_{2} v_{2}}{2} - \frac{p_{1} v_{1}}{2} + \frac{p_{1} v_{2}}{2} - \frac{p_{2} v_{1}}{2} + \frac{p_{2} v_{2}}{2}$

Let’s divide both sides by \(p_{1}\):

e = (e / p1).expand()
e

$\displaystyle - v_{1} + \frac{p_{2} v_{2}}{p_{1}} = \frac{\gamma v_{1}}{2} - \frac{\gamma v_{2}}{2} + \frac{\gamma p_{2} v_{1}}{2 p_{1}} - \frac{\gamma p_{2} v_{2}}{2 p_{1}} - \frac{v_{1}}{2} + \frac{v_{2}}{2} - \frac{p_{2} v_{1}}{2 p_{1}} + \frac{p_{2} v_{2}}{2 p_{1}}$

and then by \(v_{1}\):

e = (e / v1).expand()
e

$\displaystyle -1 + \frac{p_{2} v_{2}}{p_{1} v_{1}} = \frac{\gamma}{2} - \frac{\gamma v_{2}}{2 v_{1}} + \frac{\gamma p_{2}}{2 p_{1}} - \frac{\gamma p_{2} v_{2}}{2 p_{1} v_{1}} - \frac{1}{2} + \frac{v_{2}}{2 v_{1}} - \frac{p_{2}}{2 p_{1}} + \frac{p_{2} v_{2}}{2 p_{1} v_{1}}$

Let’s move the -1 constant from the LHS to the RHS:

e = e + 1
e

$\displaystyle \frac{p_{2} v_{2}}{p_{1} v_{1}} = \frac{\gamma}{2} - \frac{\gamma v_{2}}{2 v_{1}} + \frac{\gamma p_{2}}{2 p_{1}} - \frac{\gamma p_{2} v_{2}}{2 p_{1} v_{1}} + \frac{1}{2} + \frac{v_{2}}{2 v_{1}} - \frac{p_{2}}{2 p_{1}} + \frac{p_{2} v_{2}}{2 p_{1} v_{1}}$

Now, let’s bring the following term from the RHS to the LHS:

• \(p_{2} v_{2} / (2 p_{1} v_{1})\)

• \(\gamma p_{2} v_{2} / (2 p_{1} v_{1})\)

In order to avoid typing errors, we are going to select those expressions from the arguments that compose the RHS:

toe = table_of_expressions(e.rhs)

idx	args
0	\(\frac{1}{2}\)
1	\(\frac{\gamma}{2}\)
2	\(\frac{v_{2}}{2 v_{1}}\)
3	\(- \frac{p_{2}}{2 p_{1}}\)
4	\(\frac{\gamma p_{2}}{2 p_{1}}\)
5	\(- \frac{\gamma v_{2}}{2 v_{1}}\)
6	\(\frac{p_{2} v_{2}}{2 p_{1} v_{1}}\)
7	\(- \frac{\gamma p_{2} v_{2}}{2 p_{1} v_{1}}\)

e = e - toe[6] - toe[7]
e

$\displaystyle \frac{\gamma p_{2} v_{2}}{2 p_{1} v_{1}} + \frac{p_{2} v_{2}}{2 p_{1} v_{1}} = \frac{\gamma}{2} - \frac{\gamma v_{2}}{2 v_{1}} + \frac{\gamma p_{2}}{2 p_{1}} + \frac{1}{2} + \frac{v_{2}}{2 v_{1}} - \frac{p_{2}}{2 p_{1}}$

Let’s get rid of the constant 2 at denominators:

e = (e * 2).expand()
e

$\displaystyle \frac{\gamma p_{2} v_{2}}{p_{1} v_{1}} + \frac{p_{2} v_{2}}{p_{1} v_{1}} = \gamma - \frac{\gamma v_{2}}{v_{1}} + \frac{\gamma p_{2}}{p_{1}} + 1 + \frac{v_{2}}{v_{1}} - \frac{p_{2}}{p_{1}}$

At this point, it would be nice to factor terms sharing the same ratio. For example:

\[\begin{split}\begin{aligned} \gamma \frac{p_{2}v_{2}}{p_{1}v_{1}} + \frac{p_{2}v_{2}}{p_{1}v_{1}} &= \frac{p_{2}v_{2}}{p_{1}v_{1}} (\gamma + 1) \\ \gamma \frac{p_{2}}{p_{1}} - \frac{p_{2}}{p_{1}} &= \frac{p_{2}}{p_{1}} (\gamma - 1) \\ -\gamma \frac{v_{2}}{v_{1}} + \frac{v_{2}}{v_{1}} &= -\frac{v_{2}}{v_{1}} (\gamma - 1) \\ \end{aligned}\end{split}\]

First, let’s work on the LHS, which is an addition of two terms. If we were to explore its arguments with e.lhs.args, we would see that the index of the first term is 0, and the index of the second term is 1.

e = e.applylhs(factor)
e

$\displaystyle \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1}} = \gamma - \frac{\gamma v_{2}}{v_{1}} + \frac{\gamma p_{2}}{p_{1}} + 1 + \frac{v_{2}}{v_{1}} - \frac{p_{2}}{p_{1}}$

On the RHS we can simply collect terms:

e = e.dorhs.collect(p2 / p1).collect(v2 / v1)
e

$\displaystyle \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1}} = \gamma + 1 + \frac{v_{2} \left(1 - \gamma\right)}{v_{1}} + \frac{p_{2} \left(\gamma - 1\right)}{p_{1}}$

Note, however, that we have two opposite terms: \((1 - \gamma)\) and \((\gamma - 1)\). It would be preferable to change \((1 - \gamma) = -(\gamma - 1)\).

e = e.subs(1 - gamma, Mul(-1, gamma - 1, evaluate=False))
e

$\displaystyle \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1}} = \gamma + 1 - \frac{v_{2} \left(\gamma - 1\right)}{v_{1}} + \frac{p_{2} \left(\gamma - 1\right)}{p_{1}}$

We can now divide both sides by \(\gamma - 1\):

e = e / (gamma - 1)
e

$\displaystyle \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1} \left(\gamma - 1\right)} = \frac{\gamma + 1 - \frac{v_{2} \left(\gamma - 1\right)}{v_{1}} + \frac{p_{2} \left(\gamma - 1\right)}{p_{1}}}{\gamma - 1}$

There is a small problem to be solved on the RHS: the division was not performed term by term. Right now, the RHS is a multiplication instead of an addition:

e.rhs.func

sympy.core.mul.Mul

table_of_expressions(e.rhs)

idx	args
0	\(\frac{1}{\gamma - 1}\)
1	\(\gamma + 1 - \frac{v_{2} \left(\gamma - 1\right)}{v_{1}} + \frac{p_{2} \left(\gamma - 1\right)}{p_{1}}\)

<sympy_equation.utils.table_of_expressions object at 0x72207f313ec0>

Let’s evaluate that division term by term using divide_term_by_term, exposed by this module. This function automatically detects the deonominator:

e = e.applyrhs(divide_term_by_term)
e

$\displaystyle \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1} \left(\gamma - 1\right)} = \frac{\gamma}{\gamma - 1} + \frac{1}{\gamma - 1} - \frac{v_{2}}{v_{1}} + \frac{p_{2}}{p_{1}}$

Now we can combine the terms on the RHS where the denominator is \(\gamma - 1\).

This is the perfect opportunity to showcase another function exposed by this module: process_arguments_of_add, which is going to replace the specified terms of an addition with the results of some user-provided function. First, let’s get the indices of the terms we’d like to combine:

table_of_expressions(e.rhs, mode="args")

idx	args
0	\(\frac{1}{\gamma - 1}\)
1	\(\frac{\gamma}{\gamma - 1}\)
2	\(\frac{p_{2}}{p_{1}}\)
3	\(- \frac{v_{2}}{v_{1}}\)

<sympy_equation.utils.table_of_expressions object at 0x72207f57e180>

The terms to be replaced are located at index 0 and 1. We are going to replace them with the result of SymPy’s together:

e = e.applyrhs(lambda expr: process_arguments_of_add(expr, [0, 1], together))
e

$\displaystyle \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1} \left(\gamma - 1\right)} = \frac{\gamma + 1}{\gamma - 1} - \frac{v_{2}}{v_{1}} + \frac{p_{2}}{p_{1}}$

Let’s now bring the pressure ratio to the LHS:

e = e - p2 / p1
e

$\displaystyle - \frac{p_{2}}{p_{1}} + \frac{p_{2} v_{2} \left(\gamma + 1\right)}{p_{1} v_{1} \left(\gamma - 1\right)} = \frac{\gamma + 1}{\gamma - 1} - \frac{v_{2}}{v_{1}}$

and collect it:

e = e.dolhs.collect(p2 / p1)
e

$\displaystyle \frac{p_{2} \left(-1 + \frac{v_{2} \left(\gamma + 1\right)}{v_{1} \left(\gamma - 1\right)}\right)}{p_{1}} = \frac{\gamma + 1}{\gamma - 1} - \frac{v_{2}}{v_{1}}$

Finally, we can divide both side by the term in the parenthesys of the LHS. We can easily select it with:

toe = table_of_expressions(e.lhs)

idx	args
0	\(p_{2}\)
1	\(\frac{1}{p_{1}}\)
2	\(-1 + \frac{v_{2} \left(\gamma + 1\right)}{v_{1} \left(\gamma - 1\right)}\)

e = e / toe[2]
e

$\displaystyle \frac{p_{2}}{p_{1}} = \frac{\frac{\gamma + 1}{\gamma - 1} - \frac{v_{2}}{v_{1}}}{-1 + \frac{v_{2} \left(\gamma + 1\right)}{v_{1} \left(\gamma - 1\right)}}$