vertical shock schematic

Calculating a Shock using Newton method and Van der Waals

In this post we are going to calculate a vertical (compression) shock. Therefore, the equation system is solved using the Newton-Raphson method amongst other methods. Furthermore, we use the Van der Waals equation of state to calculate fluid properties.

Introduction to Shocks

A shock is a discontinuous change of the flow condition. This phenomena can just occur in a compressible fluid. During a shock, the static pressure rises significantly, while the velocity drops. Therefore, a vertical shock transforms a supersonic flow into a subsonic flow. Hence, this is the considered case here.

Basic equations

  • R – Specific gas constant in J/kg/K
  • p – Pressure in Pa
  • T – Temperature in K
  • v – Specific volume in m³/kg
  • rho – Density in kg/m³
  • T_c – Critical temperature in K
  • p_c – Critical pressure in Pa
  • u – Velocity in m/s
  • x – Direction in m

First of all, let’s assume the following: 1-dimensional-, stationary flow and a pipe with a constant diameter. Therefore, the governing equations can be described as follows:

Continuity equation (differential)

(1)   \begin{equation*} \frac{\partial(\rho u)}{\partial x} = 0, \end{equation*}

Continuity equation (algebraic)

(2)   \begin{equation*} u_1 \rho_1 = u_2 \rho_2, \end{equation*}

Energy conservation (differential)

(3)   \begin{equation*} \frac{\partial h}{\partial x} + u \frac{\partial u}{\partial x} = 0, \end{equation*}

Energy conservation (algebraic)

(4)   \begin{equation*} h_1 + \frac{u_1^2}{2} = h_2 + \frac{u_2^2}{2}. \end{equation*}

Furthermore, you can look up the coherence between static and stagnation state here. Let us continue with the governing equations:

Momentum conservation (differential)

(5)   \begin{equation*} u \frac{\partial (u \rho)}{\partial x} + \frac{p}{\partial x} = F_v + F_f. \end{equation*}

Finally, let us neglect the volume forces F_v e.g. gravitational acceleration and the friction forces F_f. The result is:

Momentum conservation (Algebraic)

(6)   \begin{equation*} p_1 + u_1^2 \rho_1 = p_2 + u_2^2 \rho_2. \end{equation*}

Van der Waals Equation of State

Equations of state are used to calculate fluid properties. For example, the specific volume can be calculated from pressure and temperature v=v(p, T). The Van der Waals equation is used to take the real gas behaviour into account. Hence, it is an improvement compared to the the ideal gas law. Let us have a look at the basic form of the Van der Waals equation:

(7)   \begin{equation*} \left(p+\frac{a}{v^2} \right )(v-b)-RT = 0. \end{equation*}

Here, we can calculate the variable a from:

(8)   \begin{equation*} a=\frac{27}{64}\frac{R^2T_c^2}{p_c} \end{equation*}

and b from:

(9)   \begin{equation*} b=\frac{RT_c}{8p_c}\end{equation*}

Afterwards, we rearrange it with respect to v.

(10)   \begin{equation*} pv^3 - v^2(RT+bp)+av-ab = 0. \end{equation*}

Now we can see that it is a cubic equation. Hence, we need a numerical or analytical method to solve it. Additionally, we can solve the basic equation directly for the pressure p = p(T,v):

(11)   \begin{equation*} p = \left (\frac{RT}{v-b} \right )-\left (\frac{a}{v^2} \right ) \end{equation*}

In addition, caloric state variables like enthalpy h, entropy s and inner energy u can be calculated using departure functions. The derivation is wonderfully described in this this YouTube video. Hence, the the enthalpy h for a Van der Waals gas can be calculated from:

(12)   \begin{equation*} h = h_{id} + \left (\frac{-2a}{v}+\frac{bRT}{v-b} \right). \end{equation*}

Here, h_{id} is the enthalpy of an ideal gas. Hence, we can calculate it using the specific heat capacity c_p and the temperature T:

(13)   \begin{equation*} h_{id} = \int_{T_{ref}}^{T} c_{p0}(T) d\tilde{T}. \end{equation*}

Therefore, we need to define a reference state T_{ref} to solve the integral. Remember: since the reference states can differ, a single enthalpy is not worth much, always look at enthalpy differences. Furthermore, we need temperature dependend specific heat capacity of an ideal gas c_{p0}(T). For instance, you can use a polynomial to approximate the specific heat capacity.

Do do that, we can use the Python function integrate from the scipy package. As an example, let us calculate the enthalpy difference of two thermodynamic states for water (steam). Therefore, we use

  • The critical point of water T_c = 647.1 K and p_c = 220.6e5 Pa
  • Specific gas constant R = 461.5 J/kg/K
  • Reference temperature T_{ref} = 300 K
  • Polynomial for the specific heat capacity of water c_{p0} (mind the temperature range!)
  • State 1: 500 K, 2 m³/kg
  • State 2: 600 K, 2.5 m³/kg

Hence, in Python this could look like the following.

#Author: johannes <>, 31.12.17
#License: MIT License (

import scipy.integrate as integrate

Tc = 647.1 #K
pc = 220.6e5 #Pa
R = 461.5 #J/kg/K
T_ref = 300 #K

def a():
    return (27./64.)*R**2*Tc**2./pc

def b():
    return 1./8.*R*Tc/pc/8.
def get_cp_t(T):
    #Used to approximate the ideal specific heat capacity of Water (Steam)
    #Applicable from 300 - 1000 K
    return -0.0000001989*T**3 + 0.0005871194*T**2 - 0.1303510876*T + 1770.6334862186

def idealEnthalpy(T):
    return integrate.quad(get_cp_t, T_ref, T)
def vanDerWaals_h_vt(v, T):
    hIdeal = idealEnthalpy(T)[0]
    return (-2.*a()/v+(b()*R*T)/(v-b())) + hIdeal
T1 = 500 #K
T2 = 600 #K
v1 = 2   #m^3/kg
v2 = 2.5 #m^3/kg

h1 = vanDerWaals_h_vt(v1, T1)
h2 = vanDerWaals_h_vt(v2, T2)

dh = h1-h2

print dh

As a result, we get the enthalpy difference dh = -184706.91 J/kg.

Solving a cubic Equation using the Newton-Raphson Method

First of all, we can use the Newton-Raphson method to solve nonlinear equations or equation systems numerically. This is achieved by finding the roots of a function. Hence, with the Newton-Raphson method the roots of a function can be successively approximated. We can describe it as follows:

(14)   \begin{equation*} x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}. \end{equation*}

Solving a cubic Equation in Python

Let us solve the cubic equation for the specific volume using the Newton-Raphson method in Python.

Therefore, we need to know, that a cubic equation has three roots. In the two-phase area we usually have 3 real roots. Thereby, the smallest root represents the volume of the liquid, while the largest root represents the volume of the vapor. The root in between is physically not relevant.

Furthermore, we can have one real root for gas or liquid phase and two imaginary roots. Hence the real root is of relevance.

Let us apply this using Python:

#Author: johannes <>, 31.12.17
#License: MIT License (

from scipy.misc import derivative
from scipy import optimize

Tc = 647.1 #K
pc = 220.6e5 #Pa
R = 461.5 #J/kg/K

def a():
    return (27./64.)*R**2*Tc**2./pc

def b():
    return 1./8.*R*Tc/pc/8.

def vanDerWaals_v_pt(v, p, T):
    return p*v**3.-(R*T+b()*p)*v**2.+a()*v-a()*b()

def getSpecVolume(vGuess, p, T, Nmax=100, eps=1e-6):

    it = 0
    fv = vanDerWaals_v_pt(vGuess, p, T)
    dfv = derivative(vanDerWaals_v_pt, vGuess, dx = 1e-6, args = (p, T,))
    vNew = vGuess - fv/dfv
    while it < Nmax:
        v = vNew
        fv = vanDerWaals_v_pt(v, p, T)
        dfv = derivative(vanDerWaals_v_pt, v, dx = 1e-10, args = (p, T,))
        vNew = v - fv/dfv
        if abs(v-vNew)/v < eps:
        it += 1
    return vNew

p = 1e5 # Pa
T = 800 # K
vguess = 3 #m³/kg

v = getSpecVolume(vguess, p, T)
vTest = optimize.newton(vanDerWaals_v_pt, vguess, args=(p, T, ))

As a result, we obtain a specific volume of 3.68758 m³/kg. Furthermore, the result from the own function getSpecVolume is confirmed by the newton function from the scipy optimize

Vertical (compression) shock

Furthermore, let us continue with the equations to describe shock phenomena. Therefore, let us derive the revelant equations to model the vertical compression shock. By inserting Eq. 2 into Eq. 4 we obtain:

(15)   \begin{equation*} h_2(T_2, v_2) - h_1 + 0.5 \cdot u_1^2 \left [\left (\frac{v_2}{v_1} \right )^2 - 1 \right ]. \end{equation*}

Furthermore, by inserting Eq. 2 into Eq. 6 we obtain:

(16)   \begin{equation*} p_2(T_2, v_2) - p_1 + \frac{u_1^2}{v_1} \left (\frac{v_2}{v_1} - 1\right). \end{equation*}

Therefore, we can solve this implicit coupled equation system for v_2 and T_2.

Solve the equation system using Python

The equation system eqSystem which contains the shock equations can be solved using the fsolve function from the scipy optimizepackage. In addition, this packages uses methods such as the Newton-Raphson.

#Author: johannes <>, 31.12.17
#License: MIT License (

import numpy as np
from scipy import optimize

def area(D):
    return math.pi*pow(D,2)*0.25

def eqSystem(x, args):
    TDown = x[0]
    vDown = x[1]

    pUp = args[0]
    hUp = args[1]
    uUp = args[2]
    vUp = args[3]
    f = np.zeros(2)
    f[0] = vanDerWaals_h_vt(vDown, TDown) - hUp+0.5*uUp**2.*((vDown/vUp)**2.-1.)
    f[1] = vanDerWaals_p_vt(vDown, TDown) - pUp+uUp**2./vUp*(vDown/vUp-1.)
    return f
p1 = 1e5 # Pa
T1 = 800 # K
v1guess = 3 #m³/kg
v1 = optimize.newton(vanDerWaals_v_pt, v1guess, args=(p1, T1, ))

mDot = 2 #kg/s
d = 0.1 #m
G = mDot/area(d)
h1 = vanDerWaals_h_vt(v1, T1)
u1 = G*v1

x0 = [1500, 10] #Guess values of T2, v2
args = [p1, h1, u1, v1] #Arguments (constant)

sol = optimize.fsolve(eqSystem, x0, args=args, xtol=1e-6, maxfev=200)

v2 = sol[1]
T2 = sol[0]
p2 = vanDerWaals_p_vt(v2, T2)
u2 = u1/v1*v2


In order to produce some output data, we need to specify some input data:

  • Fluid: water
  • Upstream pressure p_1 = 1e5 Pa
  • Upstream temperature T_1 = 800 K
  • Diameter d = 0.1 m
  • Variation of the velocity v_1 700 – 1200 m/s.

Using the input data, we can run the script and plot the following figures. In the first figure we see the upstream Mach number over the downstream Mach number and the pressure ratio.
Pressure and mach number correlations during a shock
In the second figure we see the upstream Mach number over the temperature and specific volume ratio.
Temperature and specific volume orrelations during a shock
We can see, that we have a supersonic upstream state (Ma > 1). Furthermore, by looking at the figures, we can state the following:

  • Subsonic state after the shock — Ma_1 > 1 \land Ma_2 < 1
  • Increase in pressure after the shock — p_1 < p_2
  • Increase in temperature after the shock — T_1 < T_2
  • Decreasing specific volume after the shock — v_1 > v_2
  • Increase in entropy — s_1 < s_2
  • .


A shock can occur in a compressible fluid. During a shock, the static pressure rises significantly, while the velocity drops. This means, a supersonic flow is transformed into a subsonic flow. The equation system to model a shock can be obtained from a simple energy and momentum balance and the continuity equation. In order to solve this non-linear equation system we use the Newton-Raphson method which included in Python packages.

Furthermore, the Van der Waals equation can be used to receive fluid properties. Therefore, we need to solve a cubic equation using the Newton-Raphson method. Furthermore, we need departure functions to transform thermal state variables into caloric state variables. Therefore, we need to integrate fluid properties from polynomials.

Finally, in a vertical shock the pressure and Temperature increases non linear, while the Mach number and specific volume decreases.

What did you think?

I’d like to hear what you think about this post.

Let me know by leaving a comment below.